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SECTION  I 


INTRODUCTION 

This  final  technical  report  documents  the  current  status  and 
accomplishments  of  the  research  performed  by  the  University  of 
Dayton  for  the  Air  Force  Office  of  Scientific  Research  (AFOSR) , 
Broad  Agency  Announcement  under  Contract  No.  F49620-88-C-0040 . 
The  research  documented  herein  was  conducted  by  the  University  of 
Dayton  Research  Institute  (UDRI)  as  the  prime  contractor,  with  the 
University  of  Pittsburgh,  Institute  for  Computational  Mathematics 
and  Applications  (ICMA),  as  its  subcontractor.  Although  this 
document  is  issued  as  the  final  report  for  the  one-year  research 
performed  during  the  reporting  period,  the  research  status  and 
accomplishments  represent  a  continuation  of  the  earlier  two-year 
effort  which  was  supported  by  the  AFOSR  Fast  Algorithm  Initiative 
under  Contract  No.  F49620-85-C-0137 . 

1 .  BACKGROUND 

The  overall  theme  of  the  UDRI-ICMA  joint  research  program  is 
the  efficient  computation  of  the  development  of  the  free  turbulent 
round  jet  by  the  time-dependent  Navier-Stokes  Equations.  As  dis¬ 
cussed  in  Krishnamurthy  et  al.  (1987),  the  only  feasible  approach 
for  predicting  jet  turbulence  at  present  involves  a  combination  of 
(i)  the  direct  computation  of  the  complete  equations  on  a  coarsely 
resolved  grid  (as  dictated  by  available  computing  resources)  to 
describe  the  large-scale  motion  by  means  of  large-eddy  simulation 
(LES),  (ii)  accurate  modeling  of  subgrid-scale  (SGS)  turbulence  to 
describe  the  small-scale  motion  that  is  not  explicitly  resolved, 
and  (iii)  the  proper  coupling  of  the  SGS  turbulence  model  to  the 
LES  computations.  It  is  precisely  this  three-pronged  approach 
that  has  governed  the  conduct  of  the  present  research  program. 
Accordingly,  this  report  summarizes  the  •salient  aspects  of  the 
research  at  UDRI  and  ICMA,  addressing  the  aforementioned  items. 
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The  LES  computations  have  been  based  upon  the  ALgorithms  for 
GAs  Equations  (ALGAE)  computer  code,  which  has  been  developed  by 
ICMA.  Although  it  was  recognized  that  this  procedure  did  not 
correspond  to  a  true  LES  of  circular-jet  turbulence  in  view  of  its 
two-dimensional  formulation,  the  joint  research  program  accepted 
ALGAE  as  a  baseline  computational  procedure  that  can  be  refined 
and  optimized  for  the  LES  of  the  round  jet  through  successive 
modifications.  Thus,  ICMA  has  continued  to  be  responsible  for  the 
development  of  the  LES  numerical  algorithms  and  their  efficient 
implementation  on  vector  computers.  UDRI  has  been  responsible  for 
research  in  SGS  turbulence  modeling,  asymptotic  development  of  the 
farfield  structure,  and  integration  of  ICMA  computational  efforts. 

2.  SCOPE  OF  RESEARCH 

The  research  on  SGS  turbulence  modeling  has  investigated 
one-point  closure  models  and  focused  on  the  eddy-viscosity  model 
to  facilitate  the  incorporation  of  a  variable-viscosity  capability 
in  the  ALGAE  procedure.  The  specific  eddy-viscosity  model  that 
was  suggested  for  ICMA  application  to  LES  computation  is  based  on 
the  algebraic  mixing-length  formulation  of  Launder  et  al .  (1972). 
A  crucial  aspect  of  SGS  turbulence  is  the  asymptotic  analysis  of 
the  farfield.  Free-  and  wall-shear  turbulent  flows  asymptotically 
attain  the  so-called  fully  developed  state  when  the  flow  becomes 
self -similar .  The  recent  adverse-pressure-gradient  boundary-layer 
asymptotic  analysis  of  Bush  and  Krishnamur thy  (1987)  de.monstrated 
the  viability  of  the  mixing-length  model.  Although  the  boundary- 
layer  analysis  is  relevant  to  the  case  of  ducted- jet  flowfields, 
its  results  are  not  directly  applicable  to  the  case  of  a  free  jet. 
Therefore,  present  research  has  addressed  the  asymptotic  analysis 
of  the  fully  developed  region  of  the  round  jet  to  examine  the 
applicability  of  the  eddy-viscosity  model  to  SGS  turbulence.  The 
results  of  this  analysis,  given  in  Bush  and  Kr ishnamurthy  (1988), 
suggest  that  the  mixing-length  model  is  indeed  a  good  candidate 
for  round-jet  SGS  turbulence  These  results  provide,  in  addition, 
farfield  information,  with  which  the  conditions  imposed  in  ALGAE 
computations  on  certain  artificially  introduced  pseudo  boundaries 
must  be  consistent.  Appendix  A  documents  the  round-jet  analysis. 
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Preliminary  ALGAE  computations  at  ICMA  failed  to  reproduce 
the  analytical  solution  of  the  two-dimensional  inviscid  jet.  This 
failure  was  largely  due  to  inadequate  spatial  resolution  near  the 
jet  boundary.  Subsequent  use  of  a  nonzero  and  constant  value  of 
molecular  viscosity  did  compute  the  qualitative  features  of  the 
analytical  solution  of  the  two-dimensional  jet.  Further  computa¬ 
tional  testing  of  the  ALGAE  code  has  not  addressed  the  flowfield 
of  the  round  jet  and  the  incorporation  of  the  turbulent  viscosity 
that  is  specified  by  the  mixing-length  model.  Other  computational 
aspects  addressed  by  ICMA  research  include  the  development  of  a 
theory  dealing  with  the  construction  of  hybrid-difference  methods, 
and  the  description  of  algorithms  for  the  numerical  simulation  of 
three-dimensional  flows.  These  are  documented  in  Appendix  B. 

3.  OUTLINE  OF  REPORT 

A  brief  discussion  of  the  current  status  and  accomplishments 
of  the  research  is  presented  in  Section  II.  Section  III  lists 
the  documentation  from  the  research  sponsored  under  this  program. 
Section  IV  shows  the  research  personnel  supported  by  this  program. 
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SECTION  II 

STATUS  OF  RESEARCH 


The  following  paragraphs  summarize  the  research  progress  to 
date.  Detailed  descriptions  are  available  in  the  Appendices. 

1.  ASYMPTOTIC  FARFIELD  DEVELOPMENT 

As  emphasized  in  Krishnamurthy  et  al .  (1987),  the  initial 
SGS  turbulence  modeling  effort  has  investigated  one-point  closure 
models  and  selected  the  mixing-length  eddy-viscosity  model  for  use 
in  ICMA  computations.  A  computational  validation  of  this  model  by 
the  ALGAE  procedure  remains  to  be  carried  out.  Successful  valida¬ 
tions  by  LES  calculations  of  the  nearfield  should  serve  to  provide 
benchmark  data  for  comparison  with  more  refined  SGS  models. 

Another  criterion  for  establishing  the  suitability  of  the 
SGS  turbulence  model  is  the  degree  of  success  in  its  ability  to 
predict  the  farfield.  This  requires  the  determination  of  the 
structure  of  the  turbulent  jet  far  enough  downstream  of  the  nozzle 
exit  that  there  is  no  residual  effect  of  the  initial  conditions, 
and  self-similarity  of  the  flow  is  attained.  For  this  downstream 
region,  an  asymptotic  analysis  of  the  Reynolds  time-averaged  equa¬ 
tions  and  complementary  boundary  conditions,  with  limit-process 
expansions  developed  in  the  limit  of  large  Reynolds  number,  has 
been  completed  and  the  uniformly  valid  flow  behavior,  from  the  jet 
centerline  to  the  ambient  farfield,  has  been  determined. 

The  analysis,  shown  in  Appendix  A,  reveals  the  existence, 
far  downstream  of  the  nozzle  exit,  of  a  turbulent  core  region,  an 
irrotational  exterior  region,  and  a  distinguished  intermediate 
region.  Appropriate  independent  and  dependent  variables  for  ai* 
three  regions  are  identified  and  the  self-similar  formulations 
therein  are  obtained.  The  solutions  for  the  turbulent  normal 
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stresses  and  the  mean  pressure  in  the  core  region  are  ascertained 
through  the  consideration  of  higher-order  approximations  of  the 
boundary-value  problem,  in  conjunction  with  the  modeling  of  the 
experimental  data.  The  resulting  core-region  solutions  are  not 
uniformly  valid  at  the  outer  edge  of  this  region.  To  obtain  a 
uniformly  valid  description  of  the  jet  flow,  it  is  necessary  to 
introduce,  in  addition  to  the  downstream  core  region,  a  downstream 
exterior  region,  at  the  outer  edge  of  which  the  flow  quantities 
reach  their  ambient  values;  and  a  downstream  intermediate  region, 
in  which  there  is  a  change  from  a  core-region-like  flow  to  an 
exterior-region-like  one. 

The  farfield  aspects  of  the  jet  flow  are  not  often  discussed 
in  the  literature.  Landau  and  Lifshitz  (1959)  determine  "the  mean 
flow  in  the  jet  outside  the  turbulent  region,"  which  (roughly) 
corresponds  to  the  downstream  exterior  region.  Present  analysis 
shows  that,  to  leading  order  of  approximation,  the  exterior  region 
is  a  turbulent  region:  the  flow  is  irrotational ,  yet  there  is  a 
convection — pressure-gradient — turbulent-stress  balance  in  both 
the  axial  and  radial  momentum  equations.  Whereas  the  resulting 
mean-velocity  solutions  are  essentially  those  obtained  by  Landau 
and  Lifshitz,  the  exterior-region  stress-  and  pressure-function 
solutions  represent  new  information  on  the  farfield  flow  behavior. 

Although  the  leading-order  solutions  for  the  radial  velocity 
in  the  core  and  exterior  regions  match  directly,  the  corresponding 
solutions  for  other  flow  quantities  do  not.  An  examination  of  the 
leading-order  and  higher-order  solutions  in  these  two  regions  (but 
especially  those  for  the  core  region)  suggests  the  existence  of  a 
downstream  intermediate  region.  It  is  seen  that  in  this  region 
the  leading-order  solutions  of  all  flow  quantities  match  directly 
to  all  of  the  leading-order  core-region  solutions  (in  a  nearfield 
overlap  domain)  and  also  match  directly  to  all  of  the  leading- 
order  exterior-region  solutions  (in  a  farfield  overlap  domain). 

With  the  presentation  of  the  pertinent  core-,  exterior-,  and 
intermediate-region  solutions  and  with  the  determination  of  the 
core-region/intermediate-region  matching  and  of  the  exterior- 
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region/intermediate-region  matching,  the  uniformly  valid  picture 
of  the  structure  of  the  self-similar  turbulent  axisymmetric  jet  is 
complete.  The  analytical  predictions  of  the  distributions  of  the 
axial/radial  velocities  and  the  shear  stress  in  the  core  region 
are  compared  with  the  experimental  data  of  Wygnanski  and  Fiedler 
(1969),  with  good  results.  Available  experimental  results, 
however,  do  not  extend  to  the  exterior  region. 

The  foregoing  validation  of  the  eddy-viscosity  model  in  the 
farfield  development  of  the  jet  apart,  the  asymptotic  analysis  is 
also  of  value  in  providing  the  needed  farfield  information  for  the 
LES  computations.  The  specification  of  the  boundary  conditions 
is  a  key  issue  in  the  numerical  simulation  of  the  jet  by  the  ALGAE 
code.  Whereas  the  LES  research  addresses  the  development  of  a  jet 
discharging  into  an  unbounded  domain,  the  ALGAE-based  computation 
requires  confined  flowfields  and  artificially  introduced  pseudo 
boundaries.  The  axis  of  symmetry,  of  course,  is  a  real  boundary 
and  does  not  pose  any  difficulty.  It  is  the  farfield  boundaries 
(at  large  radial  and  axial  distances)  that  cause  problems.  The 
asymptotic  analysis  of  the  fully  developed  region  provides  a  means 
to  specify  these  boundary  conditions.  The  asymptotic  farfield 
solutions  do  apply  to  the  downstream  pseudo  boundary,  provided  it 
is  at  least  10  jet  diameters  downstream  of  the  nozzle  exit  (exper¬ 
imental  evidence  suggests  that  the  fully  developed  region  occurs 
between  10  and  50  jet  diameters).  These  solutions  also  apply  to 
the  radially  outward  boundary  (i.e.,  the  top  pseudo  boundary  where 
free-slip  wall  conditions  have  been  imposed  in  the  ALGAE  computa¬ 
tions),  but  only  at  axial  distances  exceeding  10  jet  diameters. 
Thus,  it  is  essential  that  the  ALGAE-based  predictions  must  be 
consistent  with  the  asymptotic  results  for  axial  distances  exceed¬ 
ing  10  jet  diameters. 

2.  COMPUTATIONAL  CONSIDERATIONS 

ICMA  research  during  the  reporting  period  has  considered 
three  major  computational  topics,  a  detailed  discussion  of  which 
is  presented  in  Appendix  B. 
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ASYMPTOTIC  ANALYSIS  OF  THE  FULLY  DEVELOPED  REGION  OF 
AN  INCOMPRESSIBLE,  FREE,  TURBULENT,  ROUND  JET 

By  V.  B.  BUSH 

King,  Buck  &  Associates,  Inc.,  San  Diego,  CA  92110,  USA 

AMD  L.  KRISHNAMURTHY 

University  of  Dayton  Research  Institute,  Dayton,  OH  45469,  USA 


ABSTRACT 

The  structure  of  the  farfield  turbulent  region  of  an  incompressible  free 
jet  developing  downstream  of  an  axisymmetric  nozzle  is  studied  by  means  of  the 
Reynolds  time-averaged  equations.  The  analysis  employs  the  method  of  matched 
asymptotic  expansions,  with  limit-process  expansions  developed  in  the  limit  of 
large  Reynolds  number.  The  analysis  reveals  the  existence,  far  downstream  of 
the  nozzle  exit,  of  a  turbulent  core  region,  an  irrotational  exterior  region, 
and  a  distinguished  intermediate  region.  Self-similar  formulations  are  sought 
for  all  three  regions  in  terms  of  appropriate  independent  and  dependei t 
variables.  The  stress-  and  pressure-function  solutions  for  the  exterio 
region,  unlike  the  mean-velocity  solutions,  represent  new  information  on  the 
farfield  flow  behavior.  The  analytical  results  of  the  centerline  decay  of 
the  mean  axial  velocity  and  those  of  the  radial  distributions  of  the  axial  and 
radial  mean-velocity  components  and  the  shear-  and  normal-stress  components 
are  compared  with  available  experimental  data. 


1. 


INTRODUCTION 


The  self-similar  turbulent  round  jet,  because  it  is  a  relatively  simple 
turbulent  shear  flow,  has  been  the  subject  of  extensive  study,  theoretical 
(see,  e.g.,  Abramovich  1963;  Hinze  1975;  Townsend  1976;  Schlichting  1979),  as 
well  as  experimental  (see,  e.g.,  Reichardt  1941;  Hinze  &  Van  der  Hegge  Zijnen 
1949;  Wygnanski  &  Fiedler  1969).  From  the  first  study  (Tollmien  1926)  up  to 
the  present,  the  (classical)  theoretical  approach  has  been  to  concentrate  on 
the  leading-order  approximation  of  the  boundary-value  problem  for  the  down¬ 
stream  core  region  to  determine  the  solutions  for  the  mean  velocity  and  the 
turbulent  shear  stress. 

This  paper  presents  a  theoretical  study  of  the  structure  of  a  turbulent 
incompressible,  isothermal  jet  issuing  from  an  axisymmetric  nozzle.  Attention 
is  directed  to  the  flowfield  region  far  enough  downstream  of  the  nozzle  exit 
that  there  is  no  residual  effect  of  the  initial  conditions  and  self-similarity 
is  attained.  In  particular,  by  means  of  a  higher-order  asymptotic  analysis 
of  the  Reynolds  time-averaged  equations  and  complementary  boundary  conditions, 
presented  in  §  2,  the  uniformly  valid  behavior  of  the  flow  quantities,  from 
the  jet  centerline  to  the  ambient  farfield,  is  determined  for  this  downstream 
self-similar  region  (see  figure  1). 

In  ^  3,  through  the  consideration  of  higher-order  approximations  of  the 
boundary-value  problem,  in  conjunction  with  the  modeling  of  the  experimental 
data,  the  solutions  for  the  turbulent  normal  stresses  and  the  mean  pressure  in 
the  core  region  are  also  ascertained.  This  higher-order  analysis  establishes 
that  the  resulting  core-region  solutions  are  not  uniformly  valid  at  the  outer 
edge  of  this  region.  To  obtain  a  uniformly  valid  picture  of  the  flowfield,  it 
is  necessary  to  introduce,  in  addition  to  the  downstream  core  region,  a  down¬ 
stream  exterior  region,  at  the  outer  edge  of  which  the  flow  quantities  attain 
their  ambient  values;  and  a  downstream  intermediate  region,  in  which  the  flow 
changes  from  a  core-region-like  flow  to  an  exterior-region-like  one. 

The  farfield  aspects  of  the  jet  flow  are  not  often  discussed  in  the  lit¬ 
erature.  Landau  &  Lifshitz  (1959)  determine  the  "mean  flow  in  the  jet  outside 
the  turbulent  region,"  which  (roughly)  corresponds  to  the  downstream  exterior 
region.  In  §  5,  the  appropriate  scaling  of  the  variables  indicates  that,  to 
leading  order  of  approximation,  the  exterior  region  is  a  turbulent  region:  the 
flow  is  irrotational ,  yet  there  is  a  convection — pressure-gradient — turbulent- 
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EXTERIOR  REGION: 

X5  =  dX,  fg  =  dR:  4  =  (rg/Xs)  =  (R/X); 

=  ifig!  U  =  d^Ug,  V  =  d^Vg,  Q  =  d^cog, 
P=d^Pg.T  =  d4Tg,  M=dVs.  N=d4 


INTERMEDIATE  REGION: 


=  dX,  r,^  =  d^^^R:  Q  =  (r^/X|^)  =  6-^'HR/Xy, 

^  =  AoX,^+dt^,^:  U  =  d^u,^,  V  =  Q  =  d^' 

P  =  d3p^  ,  T  =  6^\,  M  =  6\,  N  =  d\. 


CORE  REGION: 


X  =  dX,  r  =  R:  n  =  (r/x)  =  d^Hmy, 
vp  =  U  =  u,  V  =  dv,  Q  =  CO, 

P  =  dp,T  =  6t,  M  =  d^  N=di/. 

JET  CENTERLINE 


Figure  1.  Schematic  diagram  of  the  asymptotic  structure 


stress  balance  in  both  the  axial  and  the  radial  momentum  equations.  Whereas 
the  resulting  mean-velocity  solutions  are  essentially  the  ones  determined  by 
Landau  &  Lifshitz,  the  exterior-region  ress-  and  pressure-function  solutions 
represent  new  information  concerning  the  farfield  behavior  of  the  flow. 

Although  the  leading-order  core-region  and  exterior-region  solutions  for 
the  radial  velocity  match  directly,  the  corresponding  solutions  for  the  other 
flow  quantities  do  not.  An  examination  of  the  leading-order  and  higher-order 
solutions  for  these  two  regions  (but  especially  those  for  the  •'ore  region) 
suggests  the  existence  of  the  downstream  intermediate  region,  formulated  and 
analyzed  in  §  4.  The  leading-order  intermediate-region  solutions  of  all  flow 
quantities  match  directly  to  all  of  the  leading-order  core-region  solutions 
(in  a  nearfield  overlap  domain)  and  also  match  directly  to  all  of  the  leading- 
order  exterior-region  solutions  (in  a  farfield  overlap  domain) . 

With  the  presentation  of  the  pertinent  solutions  for  the  core,  exterior, 
and  intermediate  regions,  and  with  the  determination  of  the  core-region/ 
intermediate-region  matching  and  of  the  exterior-region/intermediate-region 
matching,  the  uniformly  valid  description  of  the  structure  of  the  self-similar 
turbulent  axisymmetric  jet  is  complete.  The  analytical  predictions  of  the 
core-region  distributions  of  the  axial  and  radial  velocities  and  the  shear 
stress  are  compared  in  §  6  with  the  experimental  data  of  Wygnanslci  &  Fiedler, 
with  good  results.  Available  experimental  results,  however,  do  not  extend  to 
the  exterior  region,  as  defined  in  this  paper. 
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2.  EQUATIONS  OF  MEAN  MOTION 

Consider  the  steady  flow  of  the  axisymmetric/round  fully  developed 
turbulent  jet  of  a  homogeneous,  incompressible  fluid  (P  =  const.).  Let  X  = 

(X-X^)/B  and  R  =  R/B  represent  the  axial  and  radial  coordinates,  with  X,  R  = 

0  J  J  ~ 

0  denoting  the  origin  of  the  jet,  B.,  the  initial  jet  radius,  and  X^,  the 
"origin  of  similarity."  The  mean  velocity  components  in  the  axial  and  radial 
directions,  respectively,  are  U  =  U/U^  and  V  =  V/U^,  with  ,  the  reference 

jet-exit  speed;  the  mean  pressure  is  P  =  (P-P^)/PU^,  with  P^  denoting  the 

ambient  pressure.  _ The  turbulent  shear-  and  normal-stress  components  are 

T  =  =  -(u'  v'  )/U^  and  M  =  T^^  =  -(u'  ^)/U^  and  N  =  T„„  =  -(v'  ^/U^. 

XRRX  J  AA  J  KK  J 

In  the  foregoing,  the  tilde  quantities  are  dimensional,  the  primes  denote  the 
fluctuating  quantities,  and  the  overbars  denote  time-averaging. 


The  continuity  and  momentum  equations  describing  the  mean  flow  are 


au  1  a(Rv)  _  v  =  -i-— • 

^  R  aR  ■  "  R  aR  '  ''  R  ax  ' 


(2.1) 


fn  ^ 

fi  a(RT) 

^  aM^ 

*  ^  = 

R  aR 

V. 

^  5xJ 

(2.2a) 


f„  3'' . 

fi  a(RN) 

aRj 

^  5r  = 

R  ""^R 

V. 

^  5xJ 

(2.2b) 


The  farfield  and  centerline  boundary  conditions  for  (2.1)  and  (2.2)  are 


U,  V,  P,  T,  M,  N  -♦  0  as  R  -»  (2.3a) 

V,  T  -*  0,  U  -♦  finite  as  R  -»  0.  (2.3b) 

The  downstream  analysis  considered  here  does  not  address  the  initial 
conditions . 
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3. 


THE  DOVNSTREXM  CORE  REGION 


Attention  is  directed  to  the  flow  region  far  downstream  of  the  nozzle 
exit,  characterized  by  x  =  6x  and  r  =  R,  with  x,  r  ~  0(1),  and  with  the 
stretching  parameter  6  <<  1,  such  that  (R/X)  =  6  (r/x)  ~  0(6).  (In  §  6,  from  a 
comparison  of  theory  and  experiment,  it  is  determined  that  6  ~  0(10  *■ )  . )  The 
flow  quantities  for  this  downstream  region  are,  in  turn,  scaled  as 

♦  (X,  R;  ...)  =  ?t'(x,  r;  6):  U  =  u,  V  =  6v,  (3.1a) 

P(X,  R;  ...)  =  6p(x,  r;  5),  (3.1b) 

T{X,  R;  . . . )  =  5t  (x,  r;  6) , 

M(X,  R;  ...)  =  6//(x,  r;  6),  N(X,  R;  — )  =  6i-'(x,  r;  6).  (3.1c) 


Thus,  the  differential  equations  of  mean  motion  in  this  region  are 


<5u  ,  1  9  ( rv )  _  I  dtp  1  dip 

^  *  F  ~^r-  =  =  F  5F'  "  F  SiF' 


f3u  ,  Oul  .  c  fl  <3(rr)  .  j.  d/u' 

“  57*  '  5F  *  *  57*  7-5^*  ®  57 


f  av  <5vl 

dx  dr 


+  i2.= 

dr 


1  djri^)  .  c- 
r  3r  dx 


(3.2) 

(3.3a) 

(3.3b) 


The  centerline  boundary  conditions  are,  now. 


V,  0,  u  ->  finite  as  r  -»  0. 


(3.4) 


In  turn,  (3.2)-(3.4)  can  be  combined  to  give  the  following  integral  relations: 


(3.5a) 


(3.5b) 


(3.5c) 
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As  r  -»  X,  subject  to  verification  of  the  farfield  behavior,  it  is  taken 
that  (3.5b)  becomes 


+  5  (p-/^)  ]  r  dr  =-  rtuv-T]j  =  0: 


Ot'' 

[  [u^  +  d(i)-f.i)]r  dr  =  {1/2)Z,  const. 
•^0 


(3.6) 


For  the  self-similar  formulation  of  this  region,  the  independent  and 
dependent  variables  are 


«  .  X,  »  =  [l]; 

V>{x,  r;  6)  =  K  ¥{t);  6):  u  =  j  ,v  =  -^"^|^|--F'  , 


p(x,  r, •§)=»;  n(i7;6), 


(3.7) 


(3.8a) 


(3.8b) 


T(x,  r;  8)  =  i?“  <^(J7;  8)  , 


^(x,  r;  8)  »  J(T?;  8),  Kx,  r;  8)  =  K(/7;  8).  (3.8c) 


Introduction  of  (3.7)  and  (3.8)  into  (3.2)  and  (3.3)  produces 


+  8|r7Nn-J)|j'  =  0, 

|^|j?(n-K)  j  +  +  Fj^|-  -  F' jjj'  =  n. 


(3.9a) 


(3.9b) 


The  primes  in  (3.8)  and  (3.9)  denote  differentiation  with  respect  to  )7.  The 
centerline  boundary  conditions  are 


F,  F'  ,  $  0,  as  ??  -»  0. 


(3.10) 


In  this  self-similar  formulation,  (3.6)  becomes 


x  2 

^  +  8(n-j)j)7  in  =  (i/2)z. 


(3.11) 
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The  introduction  of  the  pressure-integral  function,  A,  defined  by 


r  ' 

A  =  n  dJ?  ,  such  that  IT  =  A'  , 
■^0 


(3.12) 


leads  to  the  following  self-similar  boundary-value  problem  for  the  downstream 
core  region  of  the  turbulent  round  jet: 


|)7#  +  f  [^]|  +  (A'  -J)}j  =  0. 


|/J(A'-K)-\J  +  "  ^*]}]  ^ 


F,  F'  ,  A,  $ 


(3.13a) 


(3.13b) 


(3.14) 


.X  ^ 

[[^]  ^  S(A'-J)j77  dJ?  =  (1/2)Z, 


(3.15) 


A  more  detailed  examination  of  this  self-similar  downstream  region  is 
facilitated  with  the  introduction  of  the  following  asymptotic  representations: 

G(>7;  6)  ^  (^)  +  6G^  (77)  +  ...,  with  G  =  F,  A,  J,  K.  (3.16) 

Thus,  the  zeroth-  and  first-order  boundary-value  problems  are 


ho  fo  [^]]  =  »' 

[j7(A'  -K  )-A  1  =  0:  n  =  \'  , 

L  0  o'  oj  o  O' 

F  ,  F'  A  ,  $  -*  0,  f— 1  -+  B 

O'  O'  O'  0  L'T  J  o 

r  f-F'^2 

—  r?  d77  =  (1/2)Z; 

Jo  ^ 


as  7?  -»  0, 


(3.17a) 

(3.17b) 

(3.17c) 

(3.17d' 
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(3.18a) 


*  [r]^]  “  -["''''i-’o’]' 

[” 's-'  1  >  *o  *  fo  -  f  ;]]'['>  (^] '  ] 


r,,  f;.  a,,  ♦,  -0,  [^]  -B,  as  « 


0, 


dr,  =  -  dr,. 


(3.18b) 

(3.18c) 

(3.18d) 


To  proceed,  based  on  experiment  and  theory  (as  reported  by  Hinze  1975, 
Schlichting  1979,  and  others),  the  zeroth-order  approximation  for  the  axial- 
velocity  function  is  talcen  to  be 


(1  +  c  r. 


2  2  ,  2  ~  ^0 


(1  +  k) 


2  ' 


(3.19a) 


where  k  =  c^r,^.  Introduction  of  (3.19a)  into  (3.17d)  yields  (B^/c)  =  (3Z)^'^^. 
In  what  follows,  it  is  taken  that  *  1,  and,  in  turn,  c  =  (3Z)~^'  The 
farfield  and  centerline  behaviors  for  this  function  are 


l-2k 


* ...] 


-*  0  as  k  OD, 


-♦  (l-2k  +  ...)  -»  1  as  k  -»  0. 


(3.19b) 


(3.19c) 


The  corresponding  approximation  for  the  streamfunction  is 


^  k  . 

2  2,  "  *0  (1  +  k)  '  *0  ~  <3/2)Z: 


2  ( 1  +  c  r,  ) 


Fq  ->  Aq  (1-k  + _ )  Ag  as  k 


Fq  -*  A^  k(l-k  +  ...)  -*0ask->0. 


(3.20a) 


(3.20b) 

(3.20c) 
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From  (3.19)  and  (3.20),  the  zeroth-order  radial-velocity  function  is 


ro 


B„J7(1  -  c^T)^) 


2(1  +  c%^)^ 


1_  k^^Nl-k) 
(1+k)^ 


(3.21a) 


-»  0  as  k  -»  <B, 


(3.21b) 


t-r] 

[n  oj 


-  Ir-  k^'^^  (l-3k  +  ...)-»0ask-'0. 
2  C 


(3.21c) 


The  leading-order  approximation  for  the  shear-stress  function,  in  turn,  is 


F,  rFl 


#  =  -  —  f— 1  - - - - 

"  ^  2(1  +  c\ 


'  2c  (1  ,  k)3 


(3.22a) 


^  k"®^Nl-3k'^  +  ...)  -»  0  as  k 

O  ZC 


(3.22b) 


^  ^  k^^^  (1-  3k  +  ...)  -►  0  as  k  -»  0. 

0  2  c 


(3.22c) 


The  experimental  work  of  Wygnanski  &  Fiedler  (1969)  suggests  the 
following  approximations  for  and  : 


J  - - 

0  ,  ,  2  2,2 
1  +  C  ?? 


- -,  with  0  <  a„  <  1: 

(1  +  k) 


(3.23a) 


-»  “a^k  ^  (l-2k"  ^  + _ )  -»  0  as  k  -►  <», 


(3.23b) 


Jq  -»  -a^  (l-2k  +  . . . )  -»  -a^  as  k  -»  0; 


(3.23c) 


— - — -  =  -b„  - ^ - r,  with  0  <  b„  <  a„  <  1; 

0  ,,  22,2  o,,  o  0 

(1  +  c  ^  )  (1  +  k) 


(3.24a) 


-»  -b^k' ^  (l-2k' ‘  +  ...)  -►  0  as  k  ^ 


(3.24b) 


Kq  -  -bQ(l-2k  +  ...)  -  -b^j  as  k  -*  0. 


(3.24c) 


The  evaluation  of  a^  and  b^  from  experimental  data  is  presented  in  ^  6. 
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From  (3.17b),  it  follows  that 


.  *^0  ^1  /  2  r(  1  +  k)  1  1  1 . 

=  ri~ir-)  -  nnoj  = 


(3.25a) 


^^^(l-yk  ^  + _ )  -»Oask-»’>, 

0  4c  3 


(3.25b) 


-*  k^  ^  ^  F-ffi  (k  ')  -  l  +  2k-+...l  ->0ask-»0. 

0  2c  L  J 

Thus,  the  leading-order  approximation  for  the  pressure  function  is 


(3.25c) 


(3+k)  1 

LI  k  J 

(1  +  k)  ^-1 

(3.26a) 


£.k-^[i  -  I®  k-‘  .  ...] 
Ho  -*  2^  ^ndc"  ‘ )  -  3  +  6k  +  ...j 


-*  0  as  k  -»  X', 


-»  ’c  as  k  -»  0. 


(3.26b) 


(3.26c) 


The  logarithmic  blow-up  of  the  "models"  for  A^  and  FI^  as  k  -►  0  is  noted.  No 
further  consideration  is  given  to  the  centerline  blow-up,  as  the  emphasis  of 
this  paper  is  on  the  effect  of  the  "models"  for  the  pressure  and  normal 
stresses  on  the  farfield  behavior. 

Now,  the  first-order  approximation  for  the  axial-velocity  function  is 
taken  to  be 


with  =-  (3/4)(2a^-b,)  <0 

\  A  ^  A  / 


(3.27a) 


1 1  ^  -  1 

—  B^k  (1  -  2k  +...)-*  0  as  k  -»  X, 


(3.27b) 


-♦  B^  (1  -  2k  +  ...)  -►  Bj  as  k  -*  0. 


(3.27c) 
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The  value  of  is  determined  from  the  first-order  momentum  integral  relation, 
(3.18d),  with  =  1.  The  corresponding  streamf unction  is 


^  ^  TTTkT~'  ^  = 


-  (9/8)  (2a^-b^  )Z  <  0: 


{3.28a) 


F.  -  A^(l-lc-‘  .  ...)  as  k  X., 


-»  A^k(l-k  +  ...)  -*  0  as  k  -»  0. 


(3.28b) 

(3.28c) 


The  first-order  radial-velocity  function  is 


--  f; 


(l  +  k)" 


(3.29a) 


_ i_ 

2  c 


k“‘  "  (1 


3k 


- 1 


..) 


0  as  k  -»  X, 


(3.29b) 


B 

—  k' 
2c 


(1  -  3k  + 


0  as  k 


(3.29c) 


The  resulting  first-order  approximation  for  the  shear-stress  function  is 


rF 

r^'  -1 

rF' 

0 

i 

Q 

t  =  - 

+ 

1 

n 

ri 

J 

J 

—  +  n{n  -j  )  I 
n  <)  o  j 


(^n  r 


c  2 


J’n 


(l  +  k) 


( 1  +  k) 


(a  - b. 


(l  +  k) 


{  3  /  4  )  (  2  a  -  b  ) 

()  I) 


(l  +  k) 


(3.30a) 


1  /  2  r 


4c 


(4a  -3b  )  -  i(42a  -2<^h  'v'  ‘  + 

0  I)  3  <•  •' 


0  as  k  (3.30b) 


1  /  2  r 


2c 


b^-en{k'  '  )  -  ^(2a,^+3b,,)  +  |-(10a^+3bjk  + 


0  as  k  0 . 

(3.30c) 


[3-7] 


Here,  the  approximations  for  and  are  taken  to  be 


J  =  -a  — - - r,  with  a  =  const,  (to  be  specified):  (3.31a) 

'  '  (1+k)^  ' 

-a^^k’  ^  (l-2k'  +  ...)  -  Oask-*^,  (3.31b) 

-»  -a^  (l-2k  +  ...)  ^  -a^  as  k  -*  0;  (3.31c) 

K  =  -b  — with  b  =  const,  (to  be  specified):  (3.32a) 

^  ^  (1+k)  ‘ 

-*  -b^  k'  ^  (l-2k'  ^  +  ...)  -  Oask^x,  (3.32b) 

-*  -b^  (l-2k  +...)-»  -b^  as  k  ^  0.  (3.32c) 


From  (3.18b),  it  is  determined  that,  subject  to  the  pertinent  boundary 
conditions , 


(3.33) 


"o  (l-6k  +  k^)  .  r._r(l  +  k)l  (3  +  !^.) 

5 — 


n  =  \' 

1  1 


(3.34a) 


n  -  ^k  ‘  (1  -9k‘  ‘  +  . . .: 

1  4 


— k  ^  ( 1  -•  +  . . . )  -*  0  as  k  ->  X , 


[3.34b) 


A  b  _ 

-»-^(l-9k+  ...)  +y^-[-£n(k  *)  -3  +  6k+  ...]  ->xask->0.  (3.34c' 


Higher-order  solutions  for  the  core-region  flow  quantities  are  not 
considered  here.  It  is  noted,  nevertheless,  that  a  preliminary  study  of  the 
second-order  boundary-value  problems  indicates  that  the  farfield  behavior  ot 
the  velocity  solutions  of  this  order  is  such  that  the  momentum-integral 
relation  of  (3.11)  fails.  This  failure  stems  from  the  interaction  of  the 
core-region  quantities  with  those  of  the  intermediate  region,  analyzed  in  4. 
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From  the  preceding  developments,  it  is  now  possible  to  determine  the 
farfield  behavior  of  the  solutions  in  the  limit  of  r?  -*  »,  6  -»  0,  such  that 
9  =  ^  ^  T]  0(1) .  In  this  limit,  k  =  c^??^  =  (1/3Z)6"  ^9^  -»  ®;  and 


F  ^  Fq  +  5f^  +  . . . 


{3/2)Z  -  6|(9/2)Z^e“^  +  (9/8)  (2aQ-bQ)Z  +  ...|  +  □(6^)j, 


(3,35a) 


>  rF  >,  rF  'i 
F  ^  0  c  1 

—  “  —  +  —  +  . . . 

U  J  J  J 

6^  |9Z^e"*  +  ...|  +  Q(6)j  , 


(3.35b) 


P  J  « J  1'7  1 J 

[^(3/2)Ze“^  -6|(27/2)Z^e 


(27/2) Z^O’^  +  (9/8)  (2a^-b^)Z0~^  + 


□  (S  ) 


f  2!  +  ... 


(27/2)z'’©'^  +  (9/4)  (4a^-3b^)Z^e‘^  +  ...  I  +  0(6)  ,  (3.36a) 


(3.35c) 


J  3^  +  6J^  +  ... 


-a^  |9aQZ^e“^  +  ...  I  +  0(6) 


(3.36b) 


K  ^  +  aK^  +  ... 


-a^  |9b^z^^"*  +  ...  I  +  □(a)j; 


(3.36c) 


n  n  +  an  + 

0  1 


-»-a^  j  (27/4)bQZ^O“‘  +  (9/8)Z^9~^  + 


+  Q(a) 


(3.37) 
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4.  THE  DOWNSTREAM  INTERMEDIATE  REGION 

The  results  of  (3 . 35) - {3 .37)  indicate  that,  to  ensure  uniform  validity, 
a  region  exterior  to  the  downstream  core  region  is  required.  For  this  region, 
designated,  here,  as  the  downstream  intermediate  region,  the  appropriate 
scalings  of  the  original  independent  and  dependent  variables  are 


X,  =  §X,  r,  =  5  R; 
k  k 


(4.1) 


^  R;  ^)  =Ax  +  5  ip(x,  r;  6):  U  =  6^u,  V  =  {4.2a) 

Okkkk  k  r  k 

k 


P(X,  R;  6)  =6  p^.  ,  r^. ;  5)  , 


(4.2b) 


T(X,  R;  5)  =6  r  (x_  r_-  6), 

k  k  k 

M(X,  R;  5)  =  r^^;  6),  N(X,  R;  6)  =  rj^;  6).  (4.2c) 

These  variables  are  related  to  those  of  the  core  region  through 


£.1/2 

x^  =  x,  r^.  =  o  r; 


y>  =  A^x,  +  Sv,  ,  with  A„  =  (3/2)Z:  u  =  5  u,  .  v 

O  k  k  O  k 


^1^2  f  *0,  1 

5  -_*6v^  , 

V,  1,  ✓ 


p  =  6  p^,  T  =  6  r^,  ^  =  d  V  =  5 

Introduction  of  (4.1)  and  (4.2)  into  (2.1)  and  (2.2)  produces 

=  0-  V  1 

3x,  r,  c)r,  '  ■  ^k  r,  9r.  '  w  r.  ^x.  ' 

k  k  k  k  k  k  k 

A„  9u,  r  <^u,  5p.  ri  T,  )  <3;/ 


(4.3) 


(4.4) 


(4.5a) 


Oc.  k,c2  k,  k.»^k  1  kk.s-k  ,^cun 

- OA  -  -  +  o  u  3 +  V  '3 +  3 -  —  — 3 - +  ^3 •  (4.5b) 

3  0  dr,  r,  wax,  ^  r.  or.  Ok. 

r  k  ^  k ''  k  k k  ^  k  k  k-^ 


d  (r  V  )  dr 


[4-1] 


A  self-similar  formulation  of  this  intermediate  region  is  sought  through 
the  introduction  of  the  following  variables: 


,  0  =  (r  /x  )  ; 

k  k 

(4.6) 

f-F' ^ 

rF. 

r,;  S)  =  6):  u^  =  ?'  [^J  ^  v^.  =  -  T  ' 

{4.7a) 

Pt(x^, 

r^;  6)  =  {6;  6)  s  ?‘^\'  (©;  6), 

k  k  k 

(4.7b) 

(x^  ,  r^.  ;S)  =  < 

(x^  ,  r^;  6)  =  §)  ,  r  :  6)  =  (&:  6)  .  (4.7c) 

kWk  k  kkk  k 

The  independent  variables  of  the  intermediate  region  are  related  to  those  of 
the  core  region  by 

^  =  X,  0  -  -  d^''^{r/x)  =  (4.8) 

The  resulting  self-similar  axial-  and  radial-momentum  equations  are 


(4.9b) 


Again,  it  is  appropriate  to  introduce  asymptotic  representations  for  the 
dependent  variables,  i.e., 

G  (e;5)  a  + -  with  ,  A  f  J  .  (4.10) 

k  kO  kkkkkk 

Thus,  the  zeroth-order  equations  for  the  intermediate  region  are 


+  0 


A '  - 

kO 


=  0, 


A ' 

kO 


\ 

kO 

e 


=  0. 


{4.11a) 


(4.11b) 


[4-2] 


The  solutions  of  (4.11)  that  match  to  those  of  the  (farfield)  core  region  are 


F_  =  {jA  +  A  - 
k  0  [^4  0  1  o  J 


(4.12a) 


*.0  =  -{K^“  *  ■* 


{4.13a) 


(4.13b,c) 


"wo  =  -{"o®'"  ♦  3‘>o''o«“}- 


(4.14) 


Recall  that  A^  =  (3/2)Z,  =  -(9/8) (2a^-bQ  )Z,  ... 

for  this  intermediate  region  are  not  pursued  here. 


Higher-order  solutions 


In  the  limit  of  ^  §  -»  0,  such  that  <  =  6  0  ~  0(1),  the  variables 

of  the  intermediate  region  have  the  following  behaviors: 


F,  ^  F.  „  +  ...-»  6  (1/4) +  0(6) 

k  k  0  O 


{4.15a) 


(1/2)A^  +  0(5) 


(1/4)Aq<  +  0(6) 


&  I" 

k  kO  [_ 

\  -\o  *  ■■■  -  -"["o':' 

!,  -  K,„  -  •••  --»[%<■ 


(1/2)A^^‘‘  +  0(5) 


+  0(5) 


<  +  0(5) 


(4.15b) 


(4.15c) 


(4.16a) 


(4.16b) 


(4.16c) 


n  2s  n.  +  ...  -»  -6\nr~^  +  0(5)1. 


(4.17) 
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THE  DOWNSTREAM  EXTERIOR  REGION 


The  uniformly  valid  characterization  of  the  farfield  development  of  the 
turbulent  round  jet  is  completed  through  the  introduction  of  the  downstream 
exterior  region,  wherein  the  appropriate  variables  are 


X  = 

6X,  r  -  6R; 

(5.1) 

s 

S 

♦  (X, 

R;8)  =  V>  (x  ,  r  ;5):  U  =  6^u  ,V  =  6^v  , 

s  s  s  s  $ 

(5.2a) 

P(X, 

R;  6)  =  a'^p  (x  ,  r  ;  6)  , 

(5.2b) 

T(X,  R;  6)  =  6*t  (x  ,  r  ;  6)  , 

s  s  s 

M(X,  R;  6)  =  8'^/i  (x  ,  r  ;  5),  N(X,  R;  6)  =  6*  v  (x  ,  r  ;  6).  (5.2c) 

s  s  s  $  s  s 


In  this  region,  the  vorticity  is 


OIX,  R;  6) 


'av  au^' 

W  “  dR^ 


=  6  <y  (x 

S  i 


(5.3) 


In  terms  of  these  exterior-region  variables,  the  equations  of  motion  are 


au 

s 

.  a(r  v  ) 

.1  S  s 

=  0:  u 

$ 

^  1 

8 

dip 

8 

dx 

s 

r 

s 

ar 

s 

r 

8  8 

r 

"  a^' 
8  8 

K 

s 

J. 

au  ^ 

8 

.t.  «. 

ri 

a(r  T  ) 

8  9 

a/w  > 
8 

X  _  -  - 

8 

8 

^X 

9 

r 

^  « 

ar 

s 

r 

U 

s 

8 

av  'i 
$ 

fl 

a(r  f  ) 

8  8 

ar  . 
+  _ ll 

9x 

9 

^8  5r 

ar  ■ 

8 

8 

ax  J 

8*^ 

(5.4) 


(5.5a) 


(5.5b) 


Talcing  this  exterior  region  to  be  an  irrotational  one. 


(5.3)  becomes 


(til. 

au  > 

_ S 

ra  V 

8 

a  f 

ar 

9-^ 

Lax^ 

6x  1 

9 

r,  <5rJJ 

9 


[5-1] 


The  boundary  conditions  for  these  equations  are 


u,v,p,T  ju,p  -»0asr  -♦®. 

$S  S$8S  S 


(5.7) 


Again,  a  self-similar  formulation  is  sought  for  this  region  in  terms  of 
the  following  variables: 


i?  =  X  ,  ?:  =  (r  /X  ); 

8  S  S 


(5.8) 


ii>  (X  ,  r  ;S)  =KT  (<;  6)  : 

S  S  8  8 


u  =  ? 
8 


„  iiliUJp.  -  — l-F. 

•  I-  •  (ih") 


}]■• 


(5.9a) 


p  (X  ,  r  ;  6)  =  (K:  6)  =  «;  5), 

8  8  8  8  8 


(5.9b) 


f  (x  ,  r  ;  S)  =  15  $  (i;  ;5) , 

8  8  8  8 


aMx  ,  r  ;  6)  =  (K;  d).  p  (x  ,  r  ;  5)  =  {K;  6) 

888  8  888  8 


(5.9c) 


Note  that  <  =  (r  /x  )  =  (R/X)  =6  Os  dr). 
8  8 


In  turn,  (5.5)  and  (5.6),  subject  to  (5.7),  can  be  written  as 


F'.2- 


(5.10a) 


(5.10b) 


L  '  (i+c^)  * 


}]'-■ 


(5.11) 
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An  asymptotic  analysis  of  this  region,  in  terms  of  expansion  of  the  form 

(5.12) 

s  s 

from  a  consideration  of  (5.10)  and  (5.11),  leads  to 


G  (^;6)  ^  G  ^{K)  +  ...,  with  G  =  F  ,  A  ,  $  ,  J  ,  K  , 

s  sO  s  s  s  s  s  s 


F  ^  =  T  Fd  + 
so  2  0  [_ 


.  1 


(5.13a) 


rF'  -1 

/• 

sO 

1  ^  J 

0 

-,t<M 

II 

i+f: 

^  J 

-1/2 

t 

rF' 


sO 


sO 


F' 

K  '  sO 


(s.nb.o 

[1*':']“'''].  (5.14a) 


4  *0^ 


rF '  >, 

2  ,  r  ^ 

Fa  '  _  -  J  1 

=  _ 

sO 

=  - 

0 

<A 

0 

'  J 

4  0  J 

-  1 


A' 


sO 


sO 


sO 


•sO 

c 


F' 


sO 


4  0 


1  + 


1  /  2^ 


!+<■ 


(5.14b) 


(5.14c) 


The  results  for  the  velocity  functions  of  (5.13)  are  consistent  with  those 
reported  by  Landau  &  Lifshitz  (1959).  The  results  of  (5.14)  represent  new 
information  concerning  the  behavior  of  the  stress  functions  in  the  farfield. 
Higher-order  solutions  for  the  exterior  region  are  not  pursued  here. 

The  farfield  (K  ->  '^)  behaviors  of  these  zeroth-order  functions  are 


F  ~  ^  A  Fl  +  ^'  ^  +  ••• 

so  2  o  J 


sO 


-7^"  *  ■■•] 


0, 


r'  „  -I 

S  0  _  , 

1  .  .,-1 

— X - F 

~  K  f. 

1  +  +  .  .  . 

f.  sO 

V.  J 

2  0 

0  > 

0, 


s  o  4  ^ 


1  +  .. 


0, 


»*  > 

1  1 

2,-2  1 

-  J 

~  ^  A_  r 

1  +  .  . . 

sO  sO 

v.  J 

1  4 

0  1 

.. 

\ ' 


^  0 


-- 


0, 


-- 


(5.15a) 


[5.15b) 


(5.15c) 


(5.16a) 


(5.16b) 


(5.16c) 
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It  is  noted  that  the  behaviors  of  the  pressure  and  normal-stress  functions  are 


if  (Tt  -o  ),  (2n  -3  ) 

O  0  0  o 


=  iA^,  i.e.,  (>r^  =  0. 


(5.17) 


The  nearfield  (<  -»  0)  behaviors  of  these  functions  are 


F  „  fl  +  ^  +  ...1  ^  A^: 

sO  0  4  J  ® 


(5.18a) 


"A  . 

K  j  2  o  2  ^  ^  • 


(5.18b) 


^o)  -T‘-'  *  •■•] 


(5.18c) 


(5.19a) 


['■,0  -  ■’.o)  ~  -  T  [l  - 


-  — 

4  0 ' 


[5.19b) 


r  '  'i  r 

A'  -  -  K  „  ~  : 

s  0  C  sO  0 

V.  V. 


.2^-2  \  ,  1  ^  I  ^ 

-  k  K  l--='^  +  ...  -♦  -®. 

o  ^ 


(5.19c) 


The  behaviors  of  the  pressure  and  normal-stress  functions  are 


^  o  =  +•••)'  -J  n  +...),  K  ^  ~  ^  1  +  .  .  .  )  , 

sO  sO  0  sO  O  sO  O 


if  and  (2/i„  -  %)  = 


(5.20) 


It  is  seen  that  these  functions  (i)  satisfy  the  boundary  conditions  at 
infinity,  and  (ii)  match  to  the  farfield  behaviors  of  the  intermediate-region 
functions. 
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6.  RESULTS  AND  DISCUSSION 

The  experimental  data  for  the  centerline  velocity  as  a  function  of  axial 
distance  are  expressed,  in  the  present  notation,  as 


U 

j 


_ C _ 

(X-£  )  /  (2B . ) 
o  j 


2£ 

X  ' 


with  C  =  const. 


(determined  experimentally) , 

(6.1) 


i.e.,  the  centerline  velocity  is  inversely  proportional  to  the  axial  distance 
in  the  self-similar  downstream  zone.  Hinze  &  Van  der  Hegge  Zijnen  (1949) , 
hereinafter  denoted  as  H-VdHZ,  found  C  «  5.9,  (20  ^  0.085  for  20  <  X  <  100; 
Wygnanski  &  Fiedler  (1969),  denoted  as  W-F,  found  C  5.4,  (20  ^  0.093  for 
50  <  X  <  180.  In  the  self-similar  analysis  of  §  3,  it  is  found  that  the 
centerline  velocity  is 


with  K  =  6X,  Bq  =  1,  B^  =  -{3/i)  .  (6.2) 

From  the  W-F  data,  (6B^)  «  -  0.079  (as  is  shown  later  in  this  section).  A 
comparison  of  (6.1)  and  (6.2)  indicates  that,  if  terms  of  0(6  )  are  neglected, 

(20  tl  +  (SB^)]  and/or  6  (20"‘  [1  +  (6Bj)]: 

6  «  0.086  for  (20'‘  «  0.093,  (68^)  -  0.079.  (6.3) 

For  the  purpose  of  consistency  only,  hereafter,  W-F  is  employed  as  the  basis 
for  comparison. 

The  axial-velocity-distribution  data  are  (most  often)  presented  as 


*  U 
U 


=  u. 


t  (1  + 


with  K  =  (R/X) ,  L‘ 


The  measurements  of  H-VdHZ  give  L 


(C;  L) , 

=  const. (determined  experimentally).  (6.4) 
^  63.8;  W-F  find  ^  57.8.  When  terms  of 
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0(6^)  are  neglected,  the  core-region  analysis  of  §  3  shows  that 


U*  =  n-  ,  ,  with  7?  =  6'^  (R/X),  c  =  {3Z)‘ 

Ujk  ««  22.2 

^  (1  +  C  7?  ) 

A  comparison  of  (6.4)  and  (6.5)  yields 

L^<«(3ZS^)  ^  and/or  Z  ^  (36^L^ ) 

Z  0.78,  c  =  (3Z)“^^^  0.65,  =  (3/2)Z  1.2 

for  6  0.086,  ^  57.8. 


(6.5) 


(6.6) 


The  farfield  axial-velocity  distribution  determined  by  the  zeroth-order 

2  2  “  1 

exterior-region  analysis  of  §  5,  with  (1/2)6  A^  =  (4L  )  ,  is 


*  U  ( 4L  ^ )  ^  * 

U  =—^ - >  J - =U  (^;L), 

„2,  1/2  EXT  '  ' 


(1  +  K  ) 

with  (4L^)*‘  0.00433  for  57.8. 


(6.7) 


* 

In  figure  2a,  U.  .__(<;  L) ,  given,  by  (6.4),  is  compared  with  the  data  of  W-F. 

d  V/  K  Ct 

Not  surprisingly,  this  representation  compares  well  with  the  data.  Figure  2b 

*  w 

shows  U„^__(<;  L)  and  U„^_(<;  L)  of  (6.4)  and  (6.7),  respectively,  as  well  as 
^INTER^^*  intermediate-region  representation.  Since  the  data  of  W-F 

(and  others)  do  not  extend  to  the  exterior  region,  as  defined  in  this  paper, 
it  is  not  possible  to  make  a  farfield  comparison. 

2 

When  terms  of  0(6  )  are  neglected,  the  core-region  radial-velocity 
distribution  can  be  expressed  as 


*  V 

V  =5^- 


•  «•  L) 

2  2  ,  2  ’'cor  E  '  '  ' 


(6.8) 


2(1  +  L  C  ) 


The  zeroth-order  exterior-region  radial-velocity  distribution  is  given  by 


* 

V 


(4L^)~^  [(l-K^)^^"  ->•  1]  _  .* 

2,1/2  ^E  X  T  ■ 


<(!+<) 


(6.9) 


★ 

Figure  3a  compares  (6.8),  with  the  data  of  W-F.  Again,  the 

comparison  is  good.  The  solutions  L)  ,  (C;  L)  ,  and  (^  )  Li 
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-  0.004 


J 


are  shown  in  figure  3b. 


Since,  for  the  core  region, 


T  =  - 


(u' v'  ) 


=  -  (u'vM  =  6r  = 


the  pertinent  zeroth-order  representation  for  the  shear-stress  distribution  is 


*^1u_JU_^  - < -  .  «;L). 

ul  2(1  . 


(6.10) 


This  is  essentially  the  relation  for  the  shear-stress  distribution  employed  by 
Hinze(i975).  For  the  exterior  region. 


(u' v'  )  _ 

j 


(u'  v'  )  =6*  r  , 

$  s 


and  the  corresponding  zeroth-order  shear-stress  distribution  is 


^>^-nrvT^(4L^)-^[(i  ^  1] 


with  (4L^)"^  0.0000187  for  ^  57.8. 


(6.11) 


In  figure  4a,  *core^^'  (6.10),  is  compared  with  the  data  of  W-F,  with 

good  results.  Figure  4b  displays  the  solutions 

For  the  core  region,  the  normal-stress  components  are 


M  =  -  ^■=  -  (u'  ^)  =  -  6r^J, 

U 

j 


N  =  -  r^=  -  (v'  ^)  =6v 

U 

J 
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Figure  3b.  Core-,  intermediate-,  and  exterior-region  radial-velocity  results 


(6-7] 


where,  at  the  very  least,  the  zeroth-order  approximations /mode Is,  ,  and  , 
employed  in  the  analysis  presented  in  §  3,  require  the  specification  of  the 
centerline  values,  a^  and  ,  respectively.  The  data  of  W-F,  in  conjunction 
with  these  zeroth-order  models,  lead  to  the  normal-stress  distributions 


r - 

/  '  2  . 

) J 


/  2 


(6ao) 


1/2 


Uj 


(1  +  L^<^ 


r - 


/  2 


(6b^) 


1/2 


with  0.29,  0.25.  (6.12) 

From  (6.12),  (Sa^)  0.084,  (Sb^)  0.063:  a^  «  0.98,  b^^  0.73  for  8  0.086, 

and,  in  turn, 

(6b^  )  A!  -0.079:  -0.92  for  6  0.086, 

(SA  )  =  A  (8B  )  -0.095:  A  -1.1  for  6  0.086,  Z  0.78.  (6.13) 

10  1  1 

The  evaluation  of  (5b^  )  follows  directly  from  the  normal-stress  data/model 
comparison,  without  the  specification  of  6.  Thus,  it  is  consistent  to  employ 
(6b^)  ^  -0.079  to  evaluate  S,  as  is  done  in  (6.3). 
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ABSTRACT 


This  report  deals  with  three  major  computational  topics  associated  with  the  numerical  simulation  of 
circular-jet  flows.  The  first  is  the  development  of  a  theory  surrounding  the  construction  of  hybrid  diffemce 
methods  that  preserve  weak  but  persistent  unsteady  features  of  the  flow.  The  second  is  a  theoretical  and  numeri¬ 
cal  study  of  a  classic  plane  jet  problem  as  a  means  of  deducing  the  nature  of  far-field  boundary  conditions.  The 
third  is  a  description  of  the  algorithms  underlying  the  "dual  variable"  method  for  the  numerical  simulation  of 
three  dimensional,  incompressible  flows. 
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CHAPTER  1 

WEIGHT  SELECTION  PROCEDURES  FOR  HYBRID  DIFFERENCE  METHODS 

1.  Introduction 

To  preserve  weak  but  persistent  features  of  an  unsteady  flow  (such  as  the  shedding  of  certain  vortex  struc¬ 
tures)  during  numerical  simulations  it  is  necessary  to  use  difference  methods  with  small  numerical  dissipation. 
Unfortunately,  in  the  absence  (or  near  absence)  of  naturally  occurring  dissipative  mechanisms,  it  is  precisely  the 
numerical  viscosity  that  stabilizes  the  method.  Thus,  it  is  important  to  examine  the  prospect  of  constructing 
methods  that  are  robust  (with  respect  to  their  stability  characteristics),  but  not  uniformly  overly  dissipative. 

A  conceptually  simple  way  to  attentuate  the  dissipative  effects  of  a  difference  method  is  to  create  a  hybrid 
method  by  blending  the  given  method  with  a  less  robust  but  more  accurate  one.  The  idea  is  to  design  the 
weights  used  in  the  blending  process  so  that  in  regions  where  little  numerical  dissipation  is  needed,  the  accurate 
method  is  dominant,  whereas  in  regions  requiring  significant  numerical  dissipation  to  preserve  stability  (or  cer¬ 
tain  qualitative  features  of  the  solution  such  as  monotonicity),  the  original  scheme  prevails.  Such  self-adjusting 
methods  in  computational  fluid  dynamics  were  apparently  first  considered  by  Harten  and  Zwas  [1972],  and  a 
comprehensive  account  of  the  steps  involved  in  the  design  of  these  methods  is  contained  in  the  thesis  of  Wildo's 
[1983]. 

Hybridization  is  also  the  notion  behind  the  Flux  Corrected  Transport  (FCT)  schemes  of  Boris  and  Book 
[1973,  1976]  and  Book,  Boris,  and  Hain  [1975].  The  FCT  schemes,  originally  developed  for  one  space  dimen¬ 
sion,  were  given  a  nontrivial  multidimensional  generalization  by  Zalasek  [1979].  He  also  showed  that  they 
could  be  interpreted  in  terms  of  convex  combinations  of  flux  terms  related  to  low-order  (sffongly  dissipative) 
and  high-order  (marginally  dissipative)  difference  methods.  The  FCT  weight  selection  process  uses  a  "monotoni¬ 
city  constraint"  on  the  numerical  solution  and  has  a  particularly  simple  formulation  in  one  dimension. 
Specifically,  the  weights  are  determined  so  as  to  maximize  the  effect  of  the  high-order  method’s  flux  terms,  sub¬ 
ject  to  the  condition  that  over  any  timestep  no  extrema  are  introduced  that  would  not  also  be  present  in  the 
low-order  solution  at  the  new  time.  This  implies  that  the  total  variation  of  the  hybrid  grid  function  does  ii  t 
exceed  that  of  the  low-order  grid  function.  The  same  idea  has  been  used  to  define  tcal-variation-diminishing 
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(TVD)  schemes  (Harten  [19831.  fl984)). 

The  monotonicity  constraint  is  consistent  with  the  behavior  of  a  solution  for  the  one-dimensional  c^aiar 
convection  equation  (the  total  variation  of  such  a  solution  does  not  increase  in  time),  and  in  this  case  the  FCT 
algorithms  perform  impressively.  However,  the  situation  is  different  for  systems  of  nonlinear  conservation  laws. 
According  to  Woodward  and  Colella  [1984]: 

Then  no  such  monotonicity  constraint  is  implied  by  the  differential 
equations,  and  the  use  of  such  a  constraint  can  lead  to  difficulties. 

In  particular,  a  smooth  region  with  strong  gradients  can  be  turned 
into  a  sequence  of  discontinuous  jumps,  with  the  appearance  of  a 
staircase. 

The  unsuitability  of  a  TVD  condition  in  the  design  of  difference  methods  for  multidimensional  quasilinear  sys¬ 
tems  is  also  indicated  by  a  result  of  Rauch  [1986].  There  it  is  shown  that  unless  the  commutators  of  all  of  the 
Jacobian  matrices  appearing  in  the  system  vanish,  no  multiple  of  the  W*-*  seminorm  of  the  initial  condition  can 
bound  this  seminorm  at  a  later  dme.  Thus,  it  is  unlikely  that  a  numerical  solution  with  a  TV.  nroperty  will 
converge  to  a  solution  of  the  original  system. 

In  view  of  these  difficulties  with  the  TVD  condition,  it  is  appropriate  to  consider  other  weight  selection 
criteria.  One  such  alternative  is  based  on  the  ability  of  the  hybrid  method  to  conserve  (or  nearly  conserve)  the 
discrete  energy  of  the  numerical  solution  that  it  produces. 

In  this  chapter  we  study  such  hybrid  difference  methods  for  the  linear  convection  equation, 

u,  +  a(x)«;,  =  0 , 0  <  j:  <  i  .  t  >  0  ,  (I.l) 

subject  to  the  initial  condition, 

«(x.0)  =  \|/(^).  (I.2i 

We  assume  that  y  e  C '(-«’.“’)  and  that  a(jc)  is  continuous  and  satisfies  |i  >  a(x)  >  0  for  some  constant  m 
Since  a  hydrid  method  is  obtained  by  forming  weighted  combinations  of  the  differences  quotients  that  define  two 
consistent  methods,  our  main  concern  is  with  the  selections  of  the  weight  used  in  the  blending  process. 


2.  The  Continuous  Energy 
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Since  we  shall  base  our  weight  selection  principles  on  the  manner  in  which  the  energy  of  a  hybrid  method 
conforms  to  that  of  (1.1),  we  begin  with  an  examination  of  the  latter. 

We  first  note  that  by  using  the  method  of  characteristics,  we  may  write  the  exact  solution  of  (1.1),  (1.2)  as 


u(x,t)  =  ^(x,t)) , 


where  y  =y(x,i)  satisfies  the  equation 


If  we  define  the  energy  of  (1.1)  as 


=  t  . 


\ 


/(0^[/  u\x,t)dx]'^  , 
0 


(2.1) 


(2.2) 


then  it  is  clear  from  (2.1),  that  I{t)  is  bounded  when  ^  is.  Howevo-,  it  is  also  possible  for  J(t)  to  remain 
bounded  when  y  is  not  bounded. 


Theorem  2.1:  If  \ir\/  S  0,  then  —— 

di 


<0. 


Proof: 


I  I 

=  -2  J  a(x)y(y)\|/'(y)y,  dx  =-2f  a(yiy(y)y'(y)  dy  <0. 
0  0 


□ 

Difference  methods  for  (1.1),  (1.2)  are  frequently  analyzed  under  the  assumption  that  the  corresponding 
numerical  solution  is  periodic  in  the  discrete  space  variable.  Therefore,  it  is  of  interest  to  know  under  what  con¬ 
ditions  this  is  true  of  u  (x ,  r ). 


Theorem  (2.2):  Let  y  be  an  i-periodic  function.  Then  u  is  i-periodic  in  x ,  if  and  only  if  a  is  i-periodic. 
Proof:  If  u(x  +  i,  i)  =  u(x,  t),  then  by  (2.1)  y(x  +  i,  t)  =  y(x,  f)  -t-  u  Thus,  by  (2.2) 

j  j  j 


and  so 
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*+»  .  ’'■«  j 

f  _  -  f  _ 

J  /I  J  /J  * 


Differentiating  this  equation  with  respect  to  r ,  we  get 

1 


0  =  )-, 

Since  y,  -  -a O')  *  0,  we  have  a(y+i)  =  a(y). 
Conversely,  if  a(s+i)  =  aO),  then 


1 


flO+O  flO) 


f  ^ 

■'  a 


is  independent  of  x  and  the  above  steps  are  reversible. 


Now  we  assume  that  v  and  a  are  i-periodic.  In  general  the  energy  I(t)  is  not  conserved  even  though 
u(',  /)  is  periodic.  However,  with  regard  to  the  weighted  energy 

/(0  =  [J  dx]'^, 

■'o 

we  see  that 


dJ^  r  2 

dt  Jq  « 


(2.3) 


so  that  this  quantity  is  conserved.  Moreover,  since  J (0)  =  J (t)  >  ,  it  follows  that  in  the  case  of  periodi¬ 

city,  /(<)  is  bounded  by  [|j,  f  dx]'^. 

%  ^ 


3.  Hybrid  Methods 
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Consider  a  rectangular  mesh  with  uniform  x  and  t  spacings  h  and  x.  Let  v  be  a  mesh  function 
whose  value  at  the  mesh  point  (Jh,  mx)  is  denoted  by  Vj(m).  When  no  confusion  can  arise  we  shall  omit  writ¬ 
ing  the  dependence  on  j  and/or  tn . 

For  any  such  function  v  we  define  the  familiar  x  -directional  differences: 

=  V^±i, 

AjV  =  (vy+1  -  Vj)lh  , 

V,v  =  (v^  -  Vy_i)//l  , 

AxV  =  (v^+1  -  v^_i)/2/i  , 

=  (Vy+i  -  IVj  -I-  . 

as  well  as  analogous  differences  in  the  t -direction.  We  also  note  the  following  useful  identities  for  any  mesh 


functions  u  and  v ; 

V,(uAxv)  =  u5;^  +(VxU)V,v  ,  (3.1) 

V,(uv)  =  uV,v  +  (S~v)W^u  ,  (3.2) 

2vV,v  =  V^v^-t- /t(V,v)2  .  (3.3) 

and 

2vA,v  =  A,v^  -  x(A,  v)^  .  (3.4) 

With  this  notation  we  consider  the  following  difference  equation  for  (1.1): 

A,v  -(-a[(l -e)V,v -h0A,v]  =  O,  (3.5) 


where  a^(m)  =  Oy  =  a(Jh)  and  9  is  a  mesh  function  of  weights.  Equation  (3.5)  defines  a  consistent  explicit 
method  for  any  choice  of  9  and  reduces  to  the  "upwind"  or  "centered"  difference  method  when  Qjim)  is  respec¬ 
tively  0  or  1. 


To  study  the  energy  of  (3.5)  (and  close  the  system  of  difference  equations),  we  assume  that  v  is  Z.- 
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periodic  in  j,  i.e.,  vy^(m)  =  Vj{m).  In  view  of  Theorem  2.2  this  assumption  is  certainly  consistent  with  i- 
periodicity  in  a  and  when  h  =\JL.  Now  we  observe  that  (3.5)  may  be  rewritten  as 


A,  V  +  a(V,v  +  6/v)  =  0. 


(3.6) 


In  this  form  the  equation  clearly  reveals  its  "antidiffusive"  nature  when  0  >  0.  If  we  use  (3.1)  to 
Qh 

transform  the  term  —  5fv  and  then  multiply  the  result  by  2v/a,  we  get. 


2vA,v 

a 


+  2v  V,v  +  /iv  Vj,(0A,v)  -  h(V,0)v(V,v)  =  0. 


(3.7) 


Next  we  use  (3,4),  (3.3)  and  (3.2)  to  transform  the  first  three  terms  of  (3.7).  Upon  rearranging,  we  have 


a 


+  V,(v^  +  hQvA^v) 


(3.8) 


=  -^  (A, V  )^  +  /i5,"(0A,  V ) Vj  V  -  h  (V,  v)^  +  h (V, 0)v  (V,  v ). 

We  add  /iv(V,  v)^  to  both  sides  of  (3.8),  where  v  is  a  constant  and  observe  that  5,"(0A,v)V,v  =  (5,“0)(V,v)^. 
In  this  way  we  obtain 

Af 

- +  V,(v^+ fi0vA,v)  + /iv(V,v)^  (3.9) 


Letting 


and 


=  (A,v)2-  /,(i-v-S-0)(V,v)2+  /i(V,0)(vV,v). 


b  =  hV,0  , 


we  see  that 


Q  =  ^  (A,vf-  h(l  -  V  -S-Q)(V,v)\ 


+  -h  A0vA,  v)  +  Av(V,  v)^  =  Q  +  bhvV,  v  . 


(3.10) 


(3.11) 


(3.12) 
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If  we  define  the  weighted  energy  of  the  method  (3.S)  as 


J(m)  = 


then  under  suitable  conditions,  we  can  use  (3.12)  to  establish  the  bounded  nature  of  J  and  therefore  also  the 
boundedness  of  the  energy 


/(m)  = 


\j=o 


These  "suitable  conditions"  lead  directly  to  various  algorithms  for  the  weights  6. 


4.  Global  Weights 

We  assume  that  6  is  independent  of  the  space  index  j .  In  this  case  =  0,  and  if  we  choose  v  =  0,  then 
(3.12)  reduces  to 


ill 

a 


+  V,(v^  + /j0vA,v)  =  2  , 


(4.1) 


where 


2  =-^(A,v)2 -/,(!- 0)(V,v)2. 


If  we  multiply  (4.1)  by  hx,  sum  over  0<7  use  the  periodicity  of  v,  and  note  that 


A 

X  hi  =  j\m+\)  -  j\Qi)  , 


a 


we  see  that 


'j\m+\)  =  J\0)+  X  Qhi.  ■ 

Hence  the  weighted  energy  is  conserved  if 

Z  2/rT  =  0  . 

/. » 


(4.2) 


(4.3) 


(4.4) 


To  enforce  (4.4)  we  observe  that  by  (3.6), 


(4.5) 


Qlh  =  X[V,v  +  I  -  (1  -  e)(V,v)2  , 

where  X  s  zaJh  is  the  mesh  function  of  Courant  numbers.  Thus 

'^Qhx^iA  Q^  +  BQ  +  C)hh  , 
i 

where 

4  =^5:X(5,V. 

•  B  =  2)[X/i5^  +  V,v]V,v  ,  (4.6) 

/ 

c  =  z  . 

j 

Therefore,  (4.4)  follows  if  0  is  a  real  root  of  the  quadratic  equation. 

^2(6)sa  e^  +  se  +  c  =0.  (4.7) 

If  Che  "Courant  condition",  X  ^  1  holds,  then 

q2(0)  =  £(X-l)(V,v)2sO. 

/ 

and 

<72(1)  =  Z  X  (A,v)2  >  0  . 
j 

Hence  (4.7)  has  a  root  in  the  interval  [0,1]. 

If  we  insist  that  0  also  be  independent  of  the  time  index  m  (i.e.  0  =  constant),  then  we  cannot  in  general 
satisfy  (4.4).  However,  if  we  write  (4.5)  in  the  form 

Q/h  =  X[(l  -  |0)V,v  +  |0A,v]2  -  (l-0)(V,v)2  , 

let  A  =  4  x/h ,  and  use  the  elementary  inequality  (a+b)^  5  2(o^+6^,  then  we  do  have 

Q/h  <  A[(l  -  j0)V,v  +  |0A,v]2  -  (l-0)(V,v)2 


I 


Therefore,  using  the  periodicity  of  v ,  we  see  that 


£  Qhx  ^  <7,(0)  £  (V^vy^hh  , 


where 


^i(0)  =  A0^  +  (l-2A)e  +  2A- 1  .  (4.8) 

From  (4.3)  we  conclude  that  if  0  is  a  real  root  of  q  i(0)  =  0,  then  the  weighted  energy,  is  bounded  for  all  m . 
Note  that  if  A  5  -j,  then  ^i(O)  S  0  <  <7i(l),  and  so  there  is  a  root  in  [0,1]. 

We  summarize  our  findings  as  follows. 

Theorem  4.1  .•  If  A  =  <  -i-  and 

«  2A  -  1  +  (1-4A^’^ 


then  the  weighted  energy  J  of  the  hybrid  method  (3.5)  satisfies  J(m)  <  7(0)  for  all  m  S  0. 


Theorem  4.2:  If  X.  =  xaJh  <  1,  and  0(m)  is  the  root  of 


A(m)0^  +  B{m)  0  +  C(m)  =  0  , 


lying  in  [0,1],  where 


-B(m)  =  5^  [h\j{m)B}vj(m)  +  V,v^(m)]  V,vy(m), 
J 

Cim)  =  "£  (kj{m)  -  l)(.y^vj{m))^, 

/ 

then  the  weighted  energy  of  (3.5)  is  conserved,  i.e. 


J(,m)  =  JiO)  fw  all  m  >  0 
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5.  Local  Weights 

We  now  allow  0  to  depend  on  both  j  and  m .  To  bound  the  weighted  energy  in  this  case  we  use  a  tech¬ 
nique  of  Ladyzhenskaya  [1985]. 

We  choose  v  >  0,  multiply  (3.12)  by  hx,  and  sum  over  0^  j  ^  L,  0  <  n  <  m.  Upon  simplification,  the 
result  is 

P(m+l)  +  hv  '^(V,vfhx  =  P(0)+'^Qhx  +  h  £  fit.  (5.1) 

>.«  j,» 

where  b  and  Q  are  given  by  (3.10)  and  (3.11).  Now  assume  that  Ifi  I  <  p  and 

2  GAx  <  0  .  (5.2) 

/.  * 

Then, 

P{m+\)  +  hvYj  (y,vfhx^P(0)  -H  /ip  2  (5-3) 

y.*  /.» 

Applying  Cauchy’s  inequality  to  the  last  term  of  (5.3),  we  have 

P(m  +  l)  +  /iv  2  (^xvfhx  <  P(fj,  +  Ap(2  Phx)\'^  (y^v)^hx)'^  .  (5.4) 

But,  for  any  f ,  g  >  0,  v  >  0.  one  has  the  inequality 

fip/g  <  - 


and  it  follows  from  (5.4)  that 

Pim+l)  +  2  (Vxv)^/ix  <  J^(0)  +  2 


y. « 


I  . 


y. « 


(5.5) 


=  P{0)  + 


/ip^ 

2v 


Let 
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y(.m)  =  max  J(n) . 

Oi»Sm 


Then  the  right  side  of  (S.S)  is  majorized  by 


y(m+\)J{0)  +  ynm+l)(m+l)T  . 


Assuming  that  (m+l)x  5  fi,  we  deduce  that 


J\m+l)  +  ^1,  (V.vj^Ax  <  y(m+l)J(0)  +  y\m+l) 


hurt  I 


(5.6) 


This  inequality  in  turn  implies  that 


hU-^t\ 

(l--^)y(m+l)^J{0) 


V 

Therefore,  if  /i  5  /lo  and  t,  < - r,  then 

/loft 


y(ffi+l)S2;{0) 


(5.7) 


Now  we  fix  (  <  oo  and  divide  the  interval  [0,  t  ]  into  a  finite  number  of  sub-intervals,  each  of  whose 

y 

lengths  does  not  exceed - r.  Then  the  estimate  (5.7)  ^plies  on  each  sub-interval,  and  we  have  proven  the 

h(P? 

following  theorem. 

Theorem  5.1:  Let  fi  S  fio.  ^  ^  It  tmd  v  >  0.  If  (5.2)  holds,  then  for  any  t  <  <»,  there  is  a  number  K{t) 

such  that  for  all  m  satisfying  mx  ^  t, 

J(.m)<K(t)J{0)  . 

a 

In  view  of  (3.11)  and  (3.6),  it  is  clear  that  (5.2)  reduces  to  an  equality  if 

V  +  5/v  )  =  (1  -v-5 -e)'^(V,  V ),  (5.8) 

where  again  X  =  xo/fi .  This  is  a  system  of  nonlinear  difference  equations  for  0  which  is  to  hold  for 
)  =  0,  •  •  •  ,  L .  However,  due  to  the  L-periodicity  of  X  and  v ,  it  suffices  to  apply  it  for  j  =  1 ,  •  •  •  ,  L  and 


take  00  =  0£, . 
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We  now  examine  circumstances  under  which  (5.8)  yields  a  scheme  (3.S)  that  is  "total-variation  diminish¬ 
ing"  (TVD).  Recall  that  the  total-variation  of  a  L-periodic  mesh  function  v  may  be  defined  as 

7V(v)*  £  IV.vJ/i  , 

*=i 

and  the  difference  method 

(m +1)  =  2  (m  )  (5.9) 

k 

is  TVD  if  7y(v(m+l))  <  7V(v(m)). 

We  consider  solutions  of  (5.8)  that  satisfy 

1  -  V  >  0._i  >l-v--^.y  =  l.  •••  ,L  .  (5.10) 

Note  that  if  Xj  <  1,  then  (5.10)  admits  the  possibility  of  negative  weights. 

Theorem  5.2:  If  the  weights  0y  satisfy  (5.8)  and  (5.10),  then  (3.5)  is  TVD. 

Proof:  According  to  (5.8),  we  have 

A, V  =  -a.kXf'il  -  V  -  0y-i)'^V,v;  . 

Therefore, 

v^(ot+1)  =  (1  -  Xj'*(l  -  V  -  0y-i)'^]Vy(m)  +  -  V  -  0y_i)%>-i(m)  . 

If  (5.10)  holds,  then  Vy(m+1)  is  a  convex  combination  of  Vj{m)  and  Vj_i{m),  and  it  follows  from  a  result  of  Lax 
(See  Harten  [1984])  that  (3.5)  is  TVD. 

□ 

It  is  worth  noting  that  the  above  proof  also  shows  that 

ffu'n  [vy_i(m),vy(m)]  <  Vy(m+1)  S  max[vy_i(m),vy(ffi)]. 

Therefore,  under  the  hypothesis  of  Theorem  5.2,  if  v(m)  is  monotone  on  the  index  set  70.  '  '  '  .7i.  ificn 
v(m+l)  is  monotone  on  yVt-l,  •  •  •  ,  j\. 


1 
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The  difficulties  of  solving  system  (S.8)  may  be  avoided  if  we  consider  instead  the  uncoupled  equations, 

W,v  +  -^52  v)  =  (1  -  e)''V,v  .  (5.11) 

All  solutions  of  these  equations  are,  of  course,  contained  in  the  solutions  of  the  quadratics, 

Ae^  +  50  +  C=O,  (5.12) 

where 

A  +  /i^X<5*v)^/4, 

■B  =h  X(V,  vXS^v)  +  (V,  v)2  (5.13) 

C  =(X-lKV,v)^ 

If  we  determine  the  weights  using  (5.11),  then  we  have  no  guarantee  that  Theorem  5.1  holds.  However,  if 
we  have  a  solution  of  (5.11)  that  also  satisfies 

1  S  0  >  1  -  Y  .  (514) 

then  as  in  the  proof  of  Theorem  5.2  we  may  again  show  that  the  resulting  scheme  is  TVD  and  monotone. 

Theorem  5.3:  If  Xs  1,  then  there  is  a  solution  of  (5.11)  satisfying  (5.14).  Consequently  the  corresponding 
difference  method  is  TVD  and  monotone. 

Proof:  If  V,v  =  0,  then  we  can  take  0  =  0.  If  Vj,v  *  0,  and  we  write  (5.11)  in  the  form 

A.^l+p0)  =  (l-0)’^, 

where  p  =  hS^I(2'7^v),  then  it  is  easy  to  see  (Figure  1)  that  (5.11)  has  a  solution  O<0<lifp>-1  while  if 
p  <  -1,  then  (5.11)  has  a  solution  1  -  -  9  -  0- 

□ 
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6.  Numerical  Experiments 

Using  the  results  of  the  previous  sections,  we  can  formulate  three  algorithms  for  the  determination  of  the 
weights  that  appear  in  (3.S).  These  algorithms,  labeled  Algorithms  1,  2  and  3  in  increasing  order  of  sophistica¬ 
tion,  are  as  follows. 

Algorithm  1  (Theorem  4.1): 

1.  Choose  T  <  hl2\i. 

2.  A  =  jiT/h . 

3  0  =  [2A  -  1  +  (1^A^)'^]/2A. 

Algorithm  2  (Theorem  4.2): 

1.  Choose  t  < /i/p. 

2.  For  m  =  0,  1,  •  ■  • 

Compute  A  (m),  B  (m),  C(m)  by  (4.11). 

Solve  A  (m )0^  +  B  (ffi )0  +  C(rn )  =  0  for  0  €  [0,  1]. 

Algorithm  3  (Theorem  5.3) 

1.  Choose  T  <  h/p. 

2.  Fot  m  =0,  I,  ■  ■  ■ 

For  7  =  0 . L. 

X  =  xcUh 


Compute  A,B  ,C  by  (5.13). 

Solve  A  0^  +  B0  +  C=  Ofor0€  (1-1/X,  1]. 

These  three  algorithms  were  used  to  determine  the  weights  in  a  series  of  numerical  simulations  of  the 
square  wave  test  problems  given  in  Boris  and  Book  [1973].  For  this  problem  we  have  a(x)  =  .01  and  h  =  .01. 


If  we  take  p  =  .01,  then  the  restrictions  on  x  imposed  by  the  three  algorithms  are  all  satisfied  if  t  =  .2.  This  is 
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significantly  beneath  the  Courant  number  stability  limit  of  x  =  1  that  applies  to  the  pure  upwind  scheme  (6  a  0). 
Consequently,  as  time  progresses,  this  scheme  will  drastically  "smear"  the  shape  of  the  convected  wave.  We 
note  that  the  "exact  solution"  of  the  problem  is 

u{x,t)-  -  .Olt) ,  where 


2 

•  ± 
2 


if  n  <  X  <  n  +  .2 
otherwise  . 


n  =  0,  ±  1,±  2 . 


The  results  of  the  numerical  experiments  are  shown  in  Figure  2  after  800  time  steps  (t  =  160).  For  com¬ 
parison,  the  well  known  (and  highly  dissipative)  "upwind"  solution  is  also  included.  The  numerical  solutions  are 
shown  supenmposed  on  the  exact  solution.  It  is  obvious  that  all  three  algorithms  dramatically  reduce  the 
numerical  dissipation  of  the  upwind  scheme.  Moreover,  the  dispersive  "ripples"  of  Algorithm  1  are  substantially 
attenuated  by  Algorithm  2,  and  completely  eliminated  by  Algorithm  3. 


I 


I 
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CHAPTER  2 

THEORETICAL  AND  NUMERICAL  STUDIES  OF  A  PLANE  JET 

1.  Introduction 

The  specification  of  boundary  conditions  for  the  numerical  computations  of  fi'ee-jet  development  presents 
many  difficulties.  The  ideal  configuration  is  that  of  a  jet  discharging  into  an  unbounded  domain.  However, 
most  numerical  approaches  require  that  the  problem  be  transformed  into  a  pseudo  boundary-value  problem  on 
the  domain  of  the  computational  grid.  Thus  it  is  necessary  to  impose  far-field  conditions  on  certain  artifically 
introduced  boundaries. 

A  typical  situation  is  shown  in  Figure  1.  The  dashed  lines  represent  the  pseudo  boundaries  in  the  quarts 
plane.  Of  the  variety  of  conditions  tested,  specification  of  the  pressure  on  these  boundaries  appears  to  yield  the 
most  realistic  behavior  of  the  numerical  solution.  Appropriate  pressure  distributions  may  be  obtained  by  assum¬ 
ing  that  the  far-field  conditions  are  approximately  those  of  related  model  problems  for  which  analytical  solutions 
are  known.  For  example,  a  linear  distribution  results  from  a  Poiseuille  flow-type  assumption.  Another  possibil¬ 
ity,  which  is  the  subject  of  this  chapter,  is  to  use  the  pressure  field  corresponding  to  irrotational  flow  of  a  jet 
with  free  stream  lines.  This  is  a  classic  problem  of  a  type  first  considered  by  Helmholtz  [1868]  for  the  flow  out 
of  a  slot  (plane  jet).  Its  solution  has  been  reproduced  in  many  hydrodynamics  texts  (see,  for  example,  Milne- 
Thompson  [1955]).  The  solution  technique  uses  a  sequence  of  Schwarz-Christoffel  transformations  to  obtain  an 
implicit  representation  of  the  complex  potential  of  the  motion.  Bernoulli’s  theorem  may  then  be  employed  to 
relate  the  far-field  pressure  on  streamlines  within  the  jet  to  the  known  pressure  at  infinity.  In  the  next  section  of 
the  Chapter  we  rederive  the  solution  in  a  form  suitable  for  numerical  compulation.  Then  we  utilize  the  resulting 
exact  far-field  boundary  condition  to  provide  a  setting  of  a  problem  that  can  be  studied  numerically. 

2.  The  Plane  Jet 

Under  certain  simplifying  assumptions,  the  motion  of  a  plane  jet  issuing  from  an  aperature  in  an  infinitely 
long  wall  may  be  analytically  determined.  In  x-y  geometry,  the  assumptions  are  that  the  motion  is  steady, 
inviscid,  irrotational  and  incompressible. 


With  these  assumptions,  it  is  known  that  there  is  a  velocity  potential  <!>  and  a  stream  function  y  such  that  the 
velocity  field  {u  ,v )  is  given  by 


-20- 


V  =  -0,  =  V,  .  (2.2) 

Moreover,  since  ^  and  y  satisfy  the  Cauchy -Riemann  equations,  they  form  the  real  aiid  imaginary  parts  of  an 
analytic  function  of  a  complex  variable  z .  Thus  if  z  =  x  +  i> ,  then 

4>{xj)-n>(x.y)  =  H*,  (2.3) 

where  w  is  some  analytic  function  of  z.  The  function  w  is  called  the  complex  potential  of  the  motion.  Clearly, 
the  determination  of  w  constitutes  the  solution  of  the  problem. 

The  geometry  of  the  jet  problem  in  the  complex  z -plane  is  shown  in  Figure  2.  There  is  symmetry  about 
the  streamline  EF. 


E 


Figure  2.  z  -plane. 
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It  is  also  assumed  that  the  aperature  width  B'B  =  2a  and  that  the  uniform  speed  of  the  jet  at  C'C  is  U . 
The  curve  ABC  is  a  streamline  of  the  flow  and  in  particular  the  portion  BC  is  a  free  streamline.  This  means 
that  on  BC  (but  not  necessarily  on  AB)  the  speed  q  =  (u^+v^)’'^  is  equal  to  the  constant  U .  The  same  is  true  of 
the  curve  A'B’C .  Since  the  flow  is  inviscid,  it  is  possible  for  the  fluid  outside  the  Jet  to  remain  at  rest  (i.e.  q  is 
discontinuous  across  the  free  stream  lines).  If  o  is  the  (unknown)  coefficient  of  contraction,  i.e.  a  =  C'Clla, 
then  the  efflux  at  C  C  is  2aaU .  We  now  examine  the  consequences  of  these  assumptions  in  the  w -plane. 

Since  \|/  is  arbitrary  up  to  an  added  constant,  we  can  take  y  =  0  as  EF.  Next  we  note  that  on  C  C , 
=  V  =  -U .  Thus,  if  we  denote  the  values  of  y  at  C  and  C  by  yi  and  y2  respectively,  we  have 

(V2  -  yi)/2aa  =  -U  . 

But  the  conditions  at  E  and  the  symmeuy  imply  that  yj  = -y2-  Therefore,  on  the  streamline  AB  C  , 
y  =  yi  =  aaU  and  on  ABC,  y  =  y2  =  -caU .  Finally,  we  take  <J)=0  at  B  (and  B ').  The  condition  =  -v  =  (/ 
at  C  implies  that  0  =  -»  at  C  (and  C ).  Evidently,  0  =  =»  at  A  (and  A  ).  This  yields  the  situation  shown  in 
Figure  3. 

Next  we  observe  that 

dz  ,  dw  ._i  ^  .  ._t 

—  =  (— )  =  (0x+«  Vi)  =  i-U+lv) 

dw  dz 

=  -(u-jv)“'  =  -l/\). 


where  u  =  u-iv  is  the  complex  velocity.  If  we  write  this  in  polar  form  as 

M  =  q  , 

then 


c,  =  -U  —  =  —  =  —  e‘^ 
dw  \)  q 


(2.4) 


Note  that  u-t-tv  =  qe‘^  so  that  q  is  the  magnitude  (speed)  and  0  the  direction  of  the  velocity  vector.  Now  we 
determine  how  the  streamlines  ABC  and  ABC  transform  in  the  %  plane.  We  see  that: 


A'B'  -^\  =  Ulq  ,  0<q  <U  . 


B'C  q  =  e'*  ,  0  >  0  >  -nl  . 
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Figure  3.  w -plane. 

CB  ^  ^  >  0  2  -ji . 

BA  =  -Vlq  .  t/  2  <7  >  0  . 

Thus  we  have  the  situation  shown  in  Figure  4. 

If  we  let 

Q  =  log  ^  =  log  —  +10,  (2-5  > 

<1 

it  follows  from  Figure  3,  that  in  the  g -plane  the  streamlines  of  Figure  2  correspond  to  the  curves  shown  in  Fig 
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ising  the  thewy  of  such  transformations  (See  Milne-Thompson  [1955.  ch.  10]).  These  transformations  may  be 
constructed  to  produce  the  configuration  of  Figure  6. 


If  ft  =  nJlaaU ,  then  (2.6)  yields 
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dw  _  _1_ 

dC  ■  K 


(2.9) 


and 


(2.10) 


Combining  (2.8)  and  (2.9),  we  get 


bUdz  =  [1  + 


]dC. 


Integrating  the  equation,  we  obtain 


c 


]dC. 


(2.11) 


which  upon  evaluation  of  the  integrals  gives 

bUz  =  C.+ 1+11  +  (C^-1)’^  -  cos-*(l/0  (2.12) 

The  coefficient  of  contraction  o  is  easily  determined  from  (2.12).  Since  the  point  B  corresponds  to  C  =  1.  we 
have  bU2a  =  2+n,  which  yields  o  =  n/(n+2)  in  agreement  with  the  classical  theory. 

To  (2.8),  (2.10)  and  (2.12)  to  solve  the  jet  problem  we  may  proceed  as  follows. 


(I)  Use  (2.12)  to  trace  the  z -plane  curves  corresponding  to  the  C-plane  rays: 

^  =  p€‘“  -  7t  <  a  <  ji/2,  0  <  p  <  . 


According  to  (2.10), 


Therefore,  the  z  -plane  curves  are  the  streamlines 

y  =  (a-nJ2)lb 

lying  in  the  right  half  of  Figure  2.  In  particular  the  streamline  ABC  corresponds  to  a  =  0.  Note  that  the 
complex  inverse  cosine  is  given  by 

cos"'t  =  -i  log  (t+i(l-f^)''^)  . 
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(II)  Observe  that  if 


C  +  (C^-1)''  *  Re'P  . 


then  by  (2.8) 


Hence 


Re'P  =  u  ^  =  -  —  e‘^ 
dw  q 


q  =U/R,  e=P-(-}i, 

and  (u,  v)  the  velocity  vector  at  2(0  is  given  by, 

u  =  q  cos  0+n)  =  -q  cosP  , 


V  =  sin  (^Jr)  =  -q  sinP  . 


The  object  now  is  to  use  the  streamline  and  velocity  information  to  set  conditions  on  certain  curves  in  the 
upper  half  of  the  complex  r -plane.  Specifically,  we  wish  to  know  the  pressure  distribution  along  the  segment 
KLMN  and  the  v  velocity  compoment  on  segment  JK  in  Figure  7. 

Assume  that  streamline  RS  passes  through  P  on  segment  JK  when  ^  =  0)  =  Po^‘“  and  through  Q  on  seg¬ 
ment  LMN  when  ^  =  0  =  Pie‘“.  Using  the  previous  notation,  v^,  the  v -velocity  component  at  P  is 

U  .  a 

Vp  =  -  —  sin  Pc  . 
no 

where 

Co  +  (Co' -!)■'"  = 

Furthermore,  by  Bernoulli’s  equation,  pg ,  the  pressure  at  Q  is  given  by 

8(/^  n  1  X 


where  p,  =  pressure  at  5  (C  =  0)  and  8  =  fluid  density  are  given  constants,  and  /?  i  is  the  magnitude  of 

c.  +  (c.'-ir- 
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Figure  7,  Boundary  Curves. 

Finally,  we  note  that  the  above  discussion  applied  to  a  jet  issuing  from  a  wall.  The  solution  fi ,  9 ,  /> 
corresponding  to  a  jet  entering  the  wall  through  the  aperature  is  obtained  from  the  solution  u,v,p  of  the  previ¬ 
ous  case  by  letting  u  =-u,5  =-v,p  =  p. 

3.  Numerical  Experiments 

The  numerical  procedure  described  in  the  previous  section  was  applied  to  a  plane  jet  problem  with  an 
aperature  half  width  of  .001  (see  Figure  1).  Streamlines  of  the  anlaytic  solution  are  shown  in  Figure  8.  The 
derived  pressure  profiles  along  the  far-field  boundaries,  and  the  velocity  disuibution  at  the  mouth  of  the  jet  were 
then  used  to  set  the  boundary  conditions  for  a  series  of  numerical  simulations  of  this  problem  by  means  of  the 
computer  code  ALGAE  (Frey,  Hall  and  Porsching  [1987]).  ALGAE  produces  finite  difrerence  solutions  of  van- 
ous  two  dimensional  fluid  transport  models  under  a  wide  variety  of  geometries  and  boundary  conditions. 


Figure  8.  Analytic  Solution  Streamlines. 


Figure  9  presents  streamlines  of  a  steady  ALGAE  numerical  solution  obtained  on  a  28x35  (nonunifom) 
grid  utilizing  an  inviscid-flow  model.  Figure  10  shows  the  analogous  information  for  a  model  incorporating  a 
constant  molecular  viscosity  of  0.02.  Both  solutions  clearly  reproduce  the  qualitative  features  of  the  analytic 
solution.  The  inviscid  ALGAE  numerical  solution  shows  a  tendency  towards  recirculation  in  the  lower  right- 
hand  side  of  the  flow  region  even  though  it  is  the  inviscid  model  that  gives  rise  to  the  analytic  solution.  This  is 
due  to  a  slight  inaccuracy  in  the  resolution  of  the  jet’s  boundary  layer  where  it  meets  the  wall.  Insufficient  reso¬ 
lution  of  this  boundary  layer  leads  to  numerical  solutions  that,  as  shown  in  Figure  11,  are  completely  inaccurate. 
The  viscous  ALGAE  numerical  solution,  on  the  other  hand,  shows  no  recirculation  tendencies. 
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CHAPTER  3 

DUAL  VARIABLE  SOLUTION  OF  THREE  DIMENSIONAL.  INCOMPRESSIBLE  FLOW  PROBLEMS 

1.  Introduction 

The  dual  variable  method  (Amit,  Hall  and  Porsching  [1981],  Dougall,  Hall  and  Porsching  [1982])  is  a  dev¬ 
ice  to  reduce  the  computational  cost  of  solving  the  systems  of  (linear)  equations  that  arise  in  many  natural 
discretizations  of  fluid  transport  models.  To  date,  the  method  has  been  applied  to  two  dimensional  transient 
problems.  While  the  fundamental  ideas  of  the  method  readily  extend  to  three  dimensions,  implementations 
require  the  practical  resolution  of  certain  graph  theoretic  questions.  In  the  remaining  sections  of  this  Chapter, 
we  present  the  method  in  the  context  of  the  three  dimensional  Navier-Stokes  equations,  outline  constructions  and 
algorithms  for  its  implementation,  and  give  some  preliminary  numerical  results  and  performance  data. 

2.  The  Dual  Variable  Method 

We  describe  the  method  relative  to  a  finite  difference  discretization  of  the  well  known  three  dimensional 
Navier-Stokes  equations.  We  assume  that  the  flow  region  is  a  rectangular  parallelopiped,  three  of  whose  boun¬ 
dary  planes  form  a  portion  of  the  first  octant  of  a  Cartesian  coordinate  system.  We  do  not  give  the  finite 
difference  equations  themselves  since  they  are  obvious  generalizations  of  the  two  dimensional  MAC  equations 
presented  in  Amit,  Hall  and  Porsching  [1981]  and  Dougall,  Hall  and  Porsching  [1982]. 

To  simplify  the  presentation  in  this  section,  we  assume  that  the  boundary  conditions  are  homogeneous  and 
the  mesh  spacings  uniform.  Then  the  resulting  difference  equations  may  be  written  in  vector  form  as 

AW=0,  (2.1) 

QW  =Ai  A'^P  +  b  .  (2.2) 

Here  W  s  and  P  e  (JL  >  N)  are  vectors  of  the  unknown  velocities  and  pressures.  A  g  R^^  is  the 
incidence  matrix  of  a  directed  graph  associated  with  the  finite  difference  grid,  Q  e  R^*^  is  the  combined 
convection-diffusion  operator,  b  g  is  a  vector  of  source  terms,  and  A/  is  the  time  step. 

It  can  be  shown  that  the  rank  of  <4  is  V  -  1.  Therefore,  dim(ker(/4))  -  L  -  N  +  \.  Let  C  g 
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have  linearly  independent  columns  and  satisfy  the  orthogonality  condition  AC  =  0.  Any  C  satisfying  these  con¬ 
ditions  is  called  a  fundamental  matrix.  By  (2.1),  =  Cy  for  some  vector  of  "dual  variables"  y  €  ,  and  it 

follows  from  (2.2)  that 

C^QCy  =  C'^b  .  (2.3) 

Note  that  whereas  the  size  of  (2.1),  (2.2)  is  L  +  N,  the  size  of  (2.3)  is  L  -  N  +  1.  For  three  dimensional  prob¬ 
lems  L  =  0(iN),  so  that  for  large  N  (2.3)  is  about  1/2  the  size  of  (2.1),  (2.2). 

The  main  algorithmic  questions  associated  with  (2.3)  are  the  construction  of  C  and  the  solution  of  (2.3) 

itself. 

3.  Construction  of  C;  A  Cycle  Basis 

By  using  the  fact  that  A  is  the  incidence  matrix  of  a  digraph,  the  construction  of  C  is  straightforwored  in 
principle.  Indeed,  if  the  columns  of  C  correspond  to  any  set  of  L  -  A/  +  1  linearly  independent  cycles  of  the 
digraph,  then  C  is  a  fundamental  matrix  (Berge  and  Ghouila-Houri  [1965]).  In  the  case  of  two  dimensional 
problems  (i.e.  planar  graphs),  the  choice  of  the  elementary  cycles  leads  to  an  especially  sparse  matrix  C  (no  row 
has  more  than  two  nonzero  entries).  For  three  dimensional  problems  of  the  type  considered  here.  Ye  [1988]  has 
given  an  analogous  construction  based  on  the  four  cycle  stencils  shown  in  Figure  1.  Stencil  (a)  applys  only  on 
the  lower  horizontal  "plane"  of  nodes  in  the  graph,  stencil  (b)  on  the  interior-horizontal  planes,  while  the  stencils 
of  (c)  and  (d)  apply  respectively  on  the  vertical  "front"  and  "right"  planes.  These  cycles  are  independent  and  a 
counting  argument  shows  that  there  are  exactly  L  -  N  +  \  of  them. 

4.  Solution  of  Dual  Variable  Systems 

At  each  time  step  it  is  necessary  to  solve  the  dual  variable  system  (2.3).  In  the  two  dimensional  case  this 
was  done  directly  by  employing  a  "frontal"  method.  For  three  dimensional  problems,  the  size  of  this  system 
militates  against  such  an  approach  '  At  the  same  time,  empirical  studies  by  Hageman  [1975a,  b]  and  Mcsina 
[1988]  demonstrate  the  efficiency  of  interative  methods  in  the  solution  of  large  scale  systems  such  as  the  one 
encountered  here.  Consequently,  we  consider  the  solution  of  (2.3)  by  means  of  an  iteration. 


I  If  A  is  the  uniform  mesh  spring,  then  the  order  of  the  coefficient  matrix  C^QC  is  0  (2// )  =  0  {2Jh  ^). 


Figure  1.  Cycle  Stencils. 

Classical  iterative  methods  such  as  the  (block)  Gauss-Seidel  and  SOR  methods  are  guaranteed  to  converge 
if  the  coefficient  matrix  is  symmetric  and  posiuve  definite  (SPD).  The  same  is  true  of  projection  meJiods  such 
as  gradient  and  (preconditioned)  conjugate  gradient  methods.  Unfortunately,  the  coefficient  matrix  C^QC  in 
(2.3)  mjy  fail  to  be  SPD  by  virtue  of  a  nonuniform  mesh  and/or  convection  effects.  Thus,  we  devise  a  splitting 
of  C^QC  whose  diagonal  block  is  SPD  and  is  such  that  some  of  the  nonsymmetric  effects  of  the  original  matrix 
are  incorporated  into  this  block. 

To  describe  this  splitting,  we  let  Q  =  (^,,]  be  any  /k/  matrix  with  a  positive  diagonal, 
D  =  diag{qx\  ,  ■  ■  ,  Rjj)-  (The  matrix  Q  of  (2.3)  has  this  property).  Now  we  split  Q  as 

Q  =  D  +  E  , 

where 

f  -  [e,!  |o  i  -  j 

Let  0  >  0  and  consider  the  symmetnc  mattix, 
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0  =  O  +  e(£  +  £^)  . 


Q  will  be  strictly  diagonally  dominant  if 


qii>Yfi\qii  +  ,  i  =  1  .  •  •  •  ,J  , 


i.e.,  if 


0  <  ^  - - - -  .1=1.  •  •  •  ,J  . 

iJqi,  +  q^  I 


Thus,  if  0  satisfies  (4.1),  then  0  is  a  symmetric,  strictly  diagonally  dominant  matrix  (and  therefore  SPD). 


If  Q  is  already  symmetric  and  strictly  diagonally  dominant,  then 


qu  >  I  .  *■  =  1  .  ■ 


Therefore, 


_  qu  1  •  _  , 

2X117.71  ^  2  ■■■  • 


In  this  case,  the  choice  0  =  1/2  satisfies  (4.1)  and  in  fact  yields  Q  =  Q.  Guided  by  (4.1)  and  (4.2)  we  let 


r  =  min  r. 


Then  we  take 


1  1 
2  2 


.9r  if  r  <  - 
2 


With  this  choice  of  0  we  always  obtain  a  SPD  matrix  Q.  The  corresponding  SPD  splitting  for  (2.3)  is 


C^QC  y=C^(Q  -Q)Cy+C^b 


from  which  we  deduce  the  (outer)  iteration. 


C'QCy,  =C'^{Q  -C)Cy*-,  + 


System  (4.5)  is  of  the  form  Au  =  f  with 
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A  =  C^QC  , 


u  =  Ik  , 


and 


/  =C'^{Q  -Q)Cyk.i  +  C^b  . 


Note  that  if  0  =  0,  then  (4.5)  reduces  to  the  "Transformed  Jacobi"  method  of  Mesina  [1988].  Note  also 
that  the  factor  .9  in  (4.3)  may  be  replaced  by  any  positive  quantity  less  than  1.  Indeed,  it  can  probably  be 
replaced  by  1  in  practice.  Its  effect  is  to  guarantee  that  0  is  strictly  diagonally  dominant  (and  hence  SPD).  But 
there  are  other  conditions  that  imply  positive  definiteness.  For  example,  if  is  an  irreducible  matrix,  and  if 
r,  ^  Fj  for  some  I  <  i,  j  <J,  then  (or  Q  =  r ,  Q  is  irreducibly  diagonally  dominant  and  so  C^QC  is  SPD. 

Now  we  turn  to  the  solution  of  the  SPD  system  (4.5).  In  view  of  the  above  remarks,  it  suffices  to  con¬ 
sider  the  generic  system 

4u  =  /  ,  (4.6) 


where  4  is  an  SPD  matrix  of  order  (say)  N . 

We  apply  the  partial  preconditioned  conjugate  gradient  method  (PPCG)  to  (4.6).  The  method  can  be 

described  in  algorithmic  form  as  follows. 

Initialization: 

K  =  Number  of  PPCG  iterations  to  the  taken,  >  1. 

n  =  Number  of  preconditioned  conjugate  gradient  steps  to  be  taken  per  iteration,  n  >  1. 

Mo  =  Initial  approximation  of  u. 

P  =  SPD  splitting  or  preconditioning  matrix, 
e  =  Convergence  criterion,  e  >  0. 

1.  To  =/  -  AUq 

2.  If  rlro  <  e  then 

M  =  Mo 

Stop 
end  if 


3.  For  i  =  1  ,  ■  ,K 


Vo  =  “t-l 

For  y  =  1  ,  ■  ■  ,  n 

Solve  Pzj.i  =  r^_i  for 
If  7  =  1 ,  then 
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Pi=0 


P\  =  zo 
else 


P>  = 


7-2Z/-2 


Pi  =  Z/-1  + 
end  if 


= 


pjAPi 


vj  =  vj.i  +  ajpj 
rj  =  rj.i  -  ttjApj 

If  rjrj  <  e,  then 


Stop 
end  if 

Uk  =  V, 

4.  Print  warning  that  convergence  has  not  occurred. 

5. u=uk 

6.  Stop. 

We  remark  that  the  PPCG  algorithm  reduces  to  the  usual  preconditioned  conjugate  gradient  (PCG)  algo 
rithm  if  /f  =  1  and  n  =  N . 


Regarding  the  convergence  rate  of  the  PPCG  method,  it  can  be  shown  Luenberger  [1973]  that  if  the  eigen¬ 
values  of  the  SPD  matrix  P~^A  are 


0  <  X,/v  <  Xs-i  ^  5  ^  ^,-1  ^  ■”  ^  ^1  , 


and  e*  =  u  -  Uk,  then 


lle;kl 


lie 


0"A 


(4,7. 


where  =(X^Ax)'^.  This  proves  the  convergence  of  the  algorithm  for  any  choice  of  n.  Moreover.  il 

n  =  N ,  then  llcill,,  =  0,  i.e.  the  algorithm  converges  in  1  step  in  the  absence  of  roundoff.  This  is  a  well  kno\vn 
feature  of  conjugate  gradient  algorithms.  Since  each  iteration  of  of  the  PPCG  algorithm  incorporates  n -steps  of 
the  PCG,  we  also  have  (Luenberger  [1973]).  the  more  familiar  estimate. 


ll€*IU 


<  2 


-  yfm 


lle*_,IU 


(4K, 


Both  estimates  (4.7)  and  (4.8)  indicate  the  desirability  of  choosing  P  such  that  P  '4  has  a  narrow  spectrum 
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The  PPCG  algorithm  may  be  implemented  even  when  A  is  not  SPD.  However,  this  form  of  the  algorithm 
can  break  down  since  it  is  now  possible  for  pjApj  to  vanish.  Also  the  estimates  (4.7)  and  (4.8)  no  longer  hold. 

Each  iteration  of  the  PCG  algorithm  requires  n  system  solves  (for  zy_i)  and  n  matrix-vector  multiplica¬ 
tions  (for  Apj).  Storage  for  five  vectors  (z^_i,  pj,  Apj,  and  vj)  is  also  needed.  Clearly,  the  preconditioner  P 
should  be  chosen  so  as  to  simplify  the  solution  of  Pzy_,  =  ry_i.  For  this  reason,  the  choice 

P  =diagiau . a^v)  is  often  made.  This  gives  rise  to  the  partial  Jacobi  conjugate  gradient  (PJCG) 

method  since  it  can  be  shown  that  the  resulting  algorithm  is  a  polynomial  accelerated  version  of  the  classical 
Jacobi  method; 

=  (/  -P~^A)Xk  +  p-^f  . 

Another  possible  choice  for  P  is 

P  =(B  +  a)L)D-‘(D  -h  toZ,)^  , 

where  D  and  L  are  respectively  the  diagonal  and  lower  triangle  of  A  and  to  e  (0,2)  is  a  parameter.  This  yields 
the  partial  symmetric  succesive  over-relaxation  conjugate  gradient  (PSSORCG)  method.  Note  that  solution  of 
the  system  Pz^_i  =  amounts  to  the  solution  of  two  triangular  systems.  The  PJCG  and  PSSORCG  methods 
can  be  implemented  by  restarting  respectively  the  JCG  and  SSORCG  subroutines  of  ITPACK  software  package 
(Kincaid  et.  al.  [1981]). 

5.  Numerical  Experiments 

In  this  section  we  present  some  preliminary  results  from  a  computer  code  that  utilizes  the  algorithmic 
approach  described  in  the  first  four  sections  of  this  chapter.  This  code,  R31T,  is  under  development  at  the 
University  of  Pittsburgh  and  runs  on  the  CRAY  XMP-48  at  the  Pittsburgh  Supercomputing  Center. 

The  test  problems  involve  three  dimensional  driven  cavity  flows  as  described  in  Mansutti  et.  al.  [1987], 
The  cavity  is  a  cube  two  units  on  a  side.  No-sIip  boundary  conditions  are  applied  on  the  lid,  floor  and  the  two 
walls  that  are  transverse  to  the  motion  of  the  lid.  Depending  on  the  simulation,  either  free  or  no-slip  conditions 
hold  on  the  two  walls  that  align  with  the  lid  motion.  The  lid  translates  with  a  uniform  unit  velocity,  and  the 


Reynolds’  number  is  400. 
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Since  the  tnidplane  between  the  two  aligned  walls  is  a  plane  of  symmetry,  it  suffices  to  repalce  this  plane 
with  a  free-slip  wall  and  simulate  the  flow  in  only  one  half  of  the  cavity.  Accordingly,  a  10x10x5  uniform 
mesh  is  used  to  discretize  this  "half-problem".  Note  that  when  the  free-slip  condition  is  used  on  the  remaining 
aligned  wall,  the  problem  becomes  two-dimensional  in  the  sense  that  there  is  no  motion  in  the  transverse  direc¬ 
tion. 

Figure  2  shows  floor-to-lid  plots  of  the  driven  component  of  the  velocity  taken  at  the  vertical  centerline  of 
the  plane  of  symmetry.  The  3-D  simulation  with  a  free-slip  aligned  wall  reproduces  a  2-D  simulation  obtained 
with  the  ALGAE  code  (Frey,  Hall  and  Porsching  [1987]).  While  there  is  general  qualitative  agreement  between 
this  solution  and  that  of  Mansutti  et  al.  [1987],,  it  shows  considerably  more  dissipation  than  the  latter.  This  is 
due  to  the  well  known  numerical  diffusion  effect  introduced  by  the  use  of  upwind  differences  for  the  convection 
terms  in  R31T.  (See  Chapter  1  for  another  illustration  of  this  phenomenon.)  When  the  no-slip  condition  is  used 
on  the  aligned  wall,  the  profile  is  even  more  attenuated  in  respone  to  the  energy  dissipation  at  that  wall. 

The  10x10x5  grid  used  in  these  simulations  gives  rise  to  a  dual  variable  system  (2.3)  containing  805  unk¬ 
nowns.  This  is  not  large  enough  to  exhibit  the  advantage  of  using  an  iterative  method  for  its  solution.  In  fact,  a 
comparison  of  solution  times  shows  that  a  direct  solution  of  the  system  by  LU  decomposition  is  slightly  faster 
than  solution  by  the  PPCG  methoa  of  the  previous  section’.  This  is  not  surprising  since  other  studies  (see  for 
example,  Hageman  [1975a]  and  Mesina  [1988])  have  also  shown  that  direct  methods  are  competitive  with  and 
sometimes  faster  than  iterative  methods  for  small  problems.  Such  advantages  disappear  with  increasing  problem 
size. 


'  The  direct  solution  time  is  5  5  sec.  per  time  step  whereas  the  PPCG  method  solution  time  is  7.5  sec.  per  time  step. 
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